信号谱分析matlab实验.docx

1、信号谱分析matlab实验clearall;closeall;clc;cdF:MATLABspectralanalysis;%file path of experiment_data.mat;loadexperiment_data;figure;plot(data);N1=64;%data length;N2=256;data64=data(65:65+N1-1) data(129:129+N1-1) data(193:193+N1-1);%the three statistical samples of data using starting points of 64,128,and 192

2、 for N=64;data256=data(257:257+N2-1) data(513:513+N2-1) data(769:769+N2-1);%the three statistical samples of data using starting points of 256,512,and 768 for N=256;Ts=1;%sampling interval is 1s;我们对这个模型采样，获得1024个采样点，结合极限谱密度的图，我们可以知道我们的数据中包含三个sine波形以及附加的高斯有色噪声（colored Gaussian noise）%1.Periodogram Me

3、thods%9-6a K=8;figure(1);for time1=1:3 plot(linspace(0,1/Ts,K*N1),10*log10(abs(fft(data64(:,time1),K*N1).2/length(data64(:,time1);holdon;endholdoff;xlabel(frequency);ylabel(spectral density estimate(dB);xlim(0 (1/Ts/2);ylim(-20 20);title(9-6a periodogram,N=64,rect window,no smooth);%9-6(b) wn1=2*han

4、ning(N1);figure(2);for time1=1:3datatapering64(:,time1)=data64(:,time1).*wn1;plot(linspace(0,1/Ts,K*N1),10*log10(abs(fft(datatapering64(:,time1),K*N1).2/length(datatapering64(:,time1);holdon;endholdoff;xlabel(frequency);ylabel(spectral density estimate(dB);xlim(0 (1/Ts/2);ylim(-20 20);title(9-6b per

5、iodogram,N=64,hanning window,no smooth);%-%9-6cfigure(3);for time1=1:3Sf=(abs(fft(data64(:,time1),K*N1).2/length(data64(:,time1);for n=1:(length(Sf)-16)Sfmat(n,:)=Sf(n:n+15,:);endSfsmooth=Sfmat*ones(16,1)./16;plot(linspace(0,1/Ts,n),10*log10(Sfsmooth);holdon;endholdoff;xlabel(frequency);ylabel(spect

6、ral density estimate(dB);xlim(0 (1/Ts/2);ylim(-20 20);title(9-6c periodogram,N=64,rect window,smooth M=2);%9-6dwn1=2*hanning(N1);figure(4);for time1=1:3 datatapering64(:,time1)=data64(:,time1).*wn1; Sf=(abs(fft(datatapering64(:,time1),K*N1).2/length(datatapering64(:,time1);for n=1:(length(Sf)-16)Sfm

7、at(n,:)=Sf(n:n+15,:);endSfsmooth=Sfmat*ones(16,1)./16;plot(linspace(0,1/Ts,n),10*log10(Sfsmooth);holdon;endholdoff;xlabel(frequency);ylabel(spectral density estimate(dB);xlim(0 (1/Ts/2);ylim(-20 20);title(9-6d periodogram,N=64,hanning window,smooth M=2);从图9-6a中的周期图中可以看出，频谱假峰很多，可靠性差、分辨力低，难以从中得到正确的频谱信

8、息。其中的补零处理，只是为了减小做DFT时带来的栅栏效应，对频谱的分辨力并没有帮助。加升余弦窗处理后的频谱（图9-6b）情况好了些，这是由于矩形窗（no data tapering）的旁瓣电平起伏大，卷积运算会产生较大的失真，而升余弦窗的频谱旁瓣显著减小，由此截断得到的周期图的泄露明显减少。对比9-6c、9-6d的频率平滑处理，从频谱结果可以看到相较于没有平滑，平滑后的图形起伏明显平坦多了，假峰明显减少，可靠性得到了提升。%-%9-7afigure(5)for time2=1:3 plot(linspace(0,1/Ts,K*N2),10*log10(abs(fft(data256(:,tim

9、e2),K*N2).2/length(data256(:,time2);holdon;endholdoff;xlabel(frequency);ylabel(spectral density estimate(dB);xlim(0 (1/Ts/2);ylim(-20 20);title(9-7a periodogram,N=256,rect window,no smooth);%9-7cfigure(7);for time2=1:3Sf=abs(fft(data256(:,time2),K*N2).2/length(data256(:,time2);for n=1:(length(Sf)-16

10、)Sfmat(n,:)=Sf(n:n+15,:);endSfsmooth=Sfmat*ones(16,1)./16;plot(linspace(0,1/Ts,n),10*log10(Sfsmooth);holdon;endholdoff;xlabel(frequency);ylabel(spectral density estimate(dB);xlim(0 (1/Ts/2);ylim(-20 20);title(9-7c periodogram,N=256,rect window,smooth M=2);%-% %9-7bwn2=2*hanning(N2);figure(6);for tim

11、e2=1:3 datatapering256(:,time2)=data256(:,time2).*wn2; plot(linspace(0,1/Ts,K*N2),10*log10(abs(fft(datatapering256(:,time2),K*N2).2/length(datatapering256(:,time2);holdon;endholdoff;xlabel(frequency);ylabel(spectral density estimate(dB);xlim(0 (1/Ts/2);ylim(-20 20);title(9-7b periodogram,N=256,hanni

12、ng window,no smooth);%9-7dwn2=2*hanning(N2);figure(8);for time2=1:3 datatapering256(:,time2)=data256(:,time2).*wn2; Sf=abs(fft(datatapering256(:,time2),K*N2).2/length(datatapering256(:,time2);for n=1:(length(Sf)-16)Sfmat(n,:)=Sf(n:n+15,:);endSfsmooth=Sfmat*ones(16,1)./16;plot(linspace(0,1/Ts,n),10*l

13、og10(Sfsmooth);holdon;endholdoff;xlabel(frequency);ylabel(spectral density estimate(dB);xlim(0 (1/Ts/2);ylim(-20 20);title(9-7d periodogram,N=256,hanning window,smooth M=2)从上面四个不同处理的周期图中也可以看到，作升余弦窗处理后的周期图的泄露明显减少。频率平滑处理提高了可靠性。在处理数据点数N=64和N=256纵向比较来看，利用的数据越多，谱分析的结果可靠性越高，这符合常理。%2.sine-wave removalclear

14、;clc;clf;%to load dataload(experiment_data);% data_to_pro=data(65:128) data(129:192) data(193:256);t=64;%chosen data length!data_to_pro=data(257:512) data(513:768) data(769:1024);t=256; %chosen data length!N=8;%padding numberfor index=1:3;%to control which segment to usefreq=(0:1/(N*t-1):1); pxx=10*

15、log10(abs(fft(data_to_pro(:,index),t*N).2./t);figure(2)subplot(221)plot(freq,pxx);hold on;axis(0,0.5,-20,20);grid;title(PSD of the original data);n=(1:t);%-the first estimate-a,i=max(pxx);%estimated freq using formula 234pulse=zeros(length(pxx),1);pulse(i)=a;istore=ones(length(pxx),1);istore(i)=0;%t

16、o record pulse in order to reinsert afterwardsa=10(a/10);%estimated a is not in dBfcap=i/(t*N)*2*pi;sincap=sin(fcap.*n-(t-1)/2*fcap);%estimated sin wave in formula 233coscap=cos(fcap.*n-(t-1)/2*fcap);%estimated cos wave in formula 233pusi=-atan(sincap*data_to_pro(:,index)/(coscap*data_to_pro(:,index

17、)/pi*180;%estimate phaseacap=2*sqrt(a/t);%formula 235estimate=acap.*cos(fcap.*n+pusi/180*pi-(t-1)/2*fcap);% this is estimated cosine wave use a, f, pusi%-the first removal-removal1_add=data_to_pro(:,index)+estimate;%since the phase estimate will probably introduce a pi error,whether add or substract

18、 the estimate sine wave is unclearremoval1_sub=data_to_pro(:,index)-estimate;if(removal1_add*removal1_addremoval1_sub*removal1_sub)%correct removal will decrease the power of vector data % removal1=removal1_sub;%therefore after compare the two possible result we can decide which one to choseelse rem

19、oval1=removal1_add;endfigure(2)subplot(222)pxx=10*log10(abs(fft(removal1,t*N).2./t);plot(freq,pxx.*istore+pulse);axis(0,0.5,-20,20);hold on;gridon;title(PSD after the first removal)%-the second estimate-a,i=max(pxx);pulse(i)=a;istore(i)=0;a=10(a/10);fcap=i/(t*N)*2*pi;sincap=sin(fcap.*n-(t-1)/2*fcap)

20、;coscap=cos(fcap.*n-(t-1)/2*fcap);pusi=-atan(sincap*removal1)/(coscap*removal1)/pi*180;acap=2*sqrt(a/t);estimate=-acap.*cos(fcap.*n+pusi/180*pi-(t-1)/2*fcap);%-the second removal-removal2_add=removal1+estimate;removal2_sub=removal1-estimate;if(removal2_add*removal2_addremoval2_sub*removal2_sub) remo

21、val2=removal2_sub;else removal2=removal2_add;endfigure(2)subplot(223)pxx=10*log10(abs(fft(removal2,t*N).2./t);plot(freq,pxx.*istore+pulse);axis(0,0.5,-20,20);hold on;gridon;title(PSD after second removal)%-the third estimate-around 0.1fs-a,i=max(pxx(1:length(pxx)/8);pulse(i)=a;istore(i)=0;a=10(a/10)

22、;fcap=i/(t*N)*2*pi;sincap=sin(fcap.*n-(t-1)/2*fcap);coscap=cos(fcap.*n-(t-1)/2*fcap);pusi=-atan(sincap*removal2)/(coscap*removal2)/pi*180;acap=2*sqrt(a/t);estimate=-acap.*cos(fcap.*n+pusi/180*pi-(t-1)/2*fcap);removal3_add=removal2+estimate;removal3_sub=removal2-estimate;if(removal3_add*removal3_addr

23、emoval3_sub*removal3_sub) removal3=removal3_sub;else removal3=removal3_add;endfigure(2)subplot(224)pxx=10*log10(abs(fft(removal3,t*N).2./t);plot(freq,pxx.*istore+pulse);axis(0,0.5,-20,20);hold on;gridon;title(9-9a PSD after third removal)end可以看到用这个方法可以在很大程度上减少泄露。但是观察N=64的图，可以发现在0.2/T附近的谱线的位置还是很多变，虽然

24、相对之前的方法，我们可以很明显地看到一些分离开来的谱线，但是还是不理想。N=256的图则效果好很多， 0.2/T附近的谱线基本都是一致的，没有很大的变化，从这里也看出了样本的数目对结果还是有很大影响的。% 3.minimum-leakage spectrum estimatesM=16;f=linspace(0,0.5,1024); s=exp(1j*2*pi*(0:(M-1)*f);for time1=1:3index=N1:-1:M;for time2=1:length(index) X(time2,:)=data64(index(time2):-1:(index(time2)-M+1),

25、time1).;end temp1,temp2=size(X); R=1/temp1*(X*X);Sf(:,time1)=abs(M*1./diag(s.*(Rconj(s);endfigure;for time1=1:3 plot(linspace(0,0.5,1024),10*log10(Sf(:,time1);holdon;endylim(-30 10);xlim(0 0.5);title(9.10a minimum-leakage spectrum estimates with N=64 M=16);holdoff比较图9-10a和9-10b，可以看出相同阶数时增大数据量能提高可靠性。

26、图9-10bd可以看出滤波器阶数的增加会增大分辨力，但会使可靠性降低，高斯过程的起伏增大，0.1处的谱线峰值降低，这样我们选取滤波器阶数时需要一个折中，使频谱具有相对较好的分辨力和可靠性。对比谱分析的结果，可以看到minimum-leakage spectrum estimates的优势还是比较明显的，要好于参量方法中的AR建模方法。%4.AR MethodfunctionSf=yw_AR_zty(data,M,f)%write by for Y-W AR Method%exampledata=data64(:,1);M=16;f=linspace(0,1,1000);plot(f,Sf);f

27、=f(:);datalength=length(data);%计算自相关data=data(:);temp=xcorr(data,biased);Rx=temp(datalength:(2*datalength-1),1);%初始化b2=zeros(1,M);a=zeros(M,M);a(1,1)=-Rx(2)/Rx(1); b2(1)=(1-a(1,1)2)*Rx(1);%递推for N=2:M a(N,N)=-1*(Rx(N+1)+a(N-1,:)*fliplr(Rx(2:N) zeros(1,M-N+1)/b2(N-1);b2(N)=(1-(a(N,N)2)*b2(N-1);for p=

28、1:N-1a(N,p)=a(N-1,p)+a(N,N)*a(N-1,N-p);endendb2(M)=Rx(1)+Rx(2:(M+1)*a(M,:);%谱估计Sf=b2(M)./abs(1+a(M,:)*exp(-1j*2*pi*(1:M).*f.).2;-%9-11afigure(2);for time1=1:3M=16;f=linspace(0,1,1000);Sf=yw_AR_zty(data64(:,time1),M,f); plot(f,10*log10(Sf);holdon;endholdoff;xlabel(frequency);ylabel(spectral density estimate(dB);xlim(0 (1/Ts/2);ylim(-20 20);title(9-11a Y-W AR spectrum estimates,N=64,M=16)%9-11bfigure(3);for time2=1:3M=16;f=linspace(0,1,1000);Sf=yw_AR_zty(data256(:,time2),M,f); plot(f,10*log10(Sf);holdon;endholdoff;xlabel(frequency);ylab