www.gusucode.com > Self-adaptive Filter and Continuous Subdivision Fourier > Mymusic_rootmusic_Test.m

    %Self-adaptive Filter and Continuous Subdivision Fourier
%Transform are Listed as the Following

clear
clc
%第一种情况
load bar_1_normal_halfload;
data=bar_1_normal_halfload;
slip=0.018;
%第二种情况
%load bar_1_normal_fullload;
%data=bar_1_normal_fullload;
%slip=0.036;

fRated=50.0;


is0=data(1:10:201200,4);
is0=is0*20*sqrt(2)/4.4; 
fs=10060/10;

is0=is0-mean(is0);

is0=is0';

ts=1/fs;

wn=[40/(fs/2) 60/(fs/2)];                  

[filterb,filtera]=butter(2,wn);

is0=filter(filterb,filtera,is0);%Pay attention to adjust fs,simultaneously
is0=filter(filterb,filtera,is0);%Pay attention to adjust fs,simultaneously

is=is0(5*fs:1:20*fs-1);

isa=is';

clear is;

n100=length(isa);

is=isa(fs:1:(n100-fs-1));


M=fs;

mmm=500;

M11=M+mmm;

clear is;

is=isa;

for n=1:M11
    xn(n)=is(n+fs);
end
%% 通过不同的算法，求出xn中含有的频率值。
clear is;

is=xn';

NNN=round(length(is)*3/4);
XX=corrmtx(is,NNN,'mod');

%music算法
[sResult,wResult]=pmusic(XX,50,length(is),fs);

stem(wResult,sResult);

[vMax,mIndex]=max(sResult)


%rootmusic算法
[fResult,pResult]=rootmusic(XX,20,fs);%大于20,效果理想。20为经验值，如何确定信号的个数？

fResult(1:2:length(fResult));

pResult;

stem(fResult(1:2:length(fResult)),pResult);

fResult0=fResult(1:2:length(fResult));

ll=1;
for n=1:length(fResult0)
    if(fResult0(n)>48.0 & fResult0(n)<52.0)
        fReal(ll)=fResult0(n);
        pReal(ll)=pResult(n);
        ll=ll+1;
    end
end

stem(fReal,pReal);
axis([45 55 -inf inf]);