信号谱分析matlab实验.docx
《信号谱分析matlab实验.docx》由会员分享，可在线阅读，更多相关《信号谱分析matlab实验.docx（42页珍藏版）》请在冰豆网上搜索。
信号谱分析matlab实验
clearall;
closeall;
clc;
cdF:
\MATLAB\spectralanalysis;%filepathofexperiment_data.mat;
loadexperiment_data;
figure;
plot（data）;
N1=64;%datalength;
N2=256;
data64=[data（65:
65+N1-1）data（129:
129+N1-1）data（193:
193+N1-1）];%thethreestatisticalsamplesofdatausingstartingpointsof64,128,and192forN=64;
data256=[data（257:
257+N2-1）data（513:
513+N2-1）data（769:
769+N2-1）];%thethreestatisticalsamplesofdatausingstartingpointsof256,512,and768forN=256;
Ts=1;%samplingintervalis1s;
我们对这个模型采样，获得1024个采样点，结合极限谱密度的图，我们可以知道我们的数据中包含三个sine波形以及附加的高斯有色噪声（coloredGaussiannoise）
%1.PeriodogramMethods
%9-6a
K=8;
figure
（1）;
fortime1=1:
3plot（linspace（0,1/Ts,K*N1）,10*log10（abs（fft（data64（:
time1）,K*N1））.^2/length（data64（:
time1））））;
holdon;
end
holdoff;
xlabel（'frequency'）;
ylabel（'spectraldensityestimate（dB）'）;
xlim（[0（1/Ts/2）]）;
ylim（[-2020]）;
title（'9-6aperiodogram,N=64,rectwindow,nosmooth'）;
%9-6（b）
wn1=2*hanning（N1）;
figure
（2）;
fortime1=1:
3
datatapering64（:
time1）=data64（:
time1）.*wn1;
plot（linspace（0,1/Ts,K*N1）,10*log10（abs（fft（datatapering64（:
time1）,K*N1））.^2/length（datatapering64（:
time1））））;
holdon;
end
holdoff;
xlabel（'frequency'）;
ylabel（'spectraldensityestimate（dB）'）;
xlim（[0（1/Ts/2）]）;
ylim（[-2020]）;
title（'9-6bperiodogram,N=64,hanningwindow,nosmooth'）;
%--------------------------------------------------------------------
%9-6c
figure（3）;
fortime1=1:
3
Sf=（abs（fft（data64（:
time1）,K*N1））.^2/length（data64（:
time1）））;
forn=1:
（length（Sf）-16）
Sfmat（n,:
）=Sf（n:
n+15,:
）;
end
Sfsmooth=Sfmat*ones（16,1）./16;
plot（linspace（0,1/Ts,n）,10*log10（Sfsmooth））;
holdon;
end
holdoff;
xlabel（'frequency'）;
ylabel（'spectraldensityestimate（dB）'）;
xlim（[0（1/Ts/2）]）;
ylim（[-2020]）;
title（'9-6cperiodogram,N=64,rectwindow,smoothM=2'）;
%9-6d
wn1=2*hanning（N1）;
figure（4）;
fortime1=1:
3
datatapering64（:
time1）=data64（:
time1）.*wn1;Sf=（abs（fft（datatapering64（:
time1）,K*N1））.^2/length（datatapering64（:
time1）））;
forn=1:
（length（Sf）-16）
Sfmat（n,:
）=Sf（n:
n+15,:
）;
end
Sfsmooth=Sfmat*ones（16,1）./16;
plot（linspace（0,1/Ts,n）,10*log10（Sfsmooth））;
holdon;
end
holdoff;
xlabel（'frequency'）;
ylabel（'spectraldensityestimate（dB）'）;
xlim（[0（1/Ts/2）]）;
ylim（[-2020]）;
title（'9-6dperiodogram,N=64,hanningwindow,smoothM=2'）;
从图9-6a中的周期图中可以看出，频谱假峰很多，可靠性差、分辨力低，难以从中得到正确的频谱信息。
其中的补零处理，只是为了减小做DFT时带来的栅栏效应，对频谱的分辨力并没有帮助。
加升余弦窗处理后的频谱（图9-6b）情况好了些，这是由于矩形窗（nodatatapering）的旁瓣电平起伏大，卷积运算会产生较大的失真，而升余弦窗的频谱旁瓣显著减小，由此截断得到的周期图的泄露明显减少。
对比9-6c、9-6d的频率平滑处理，从频谱结果可以看到相较于没有平滑，平滑后的图形起伏明显平坦多了，假峰明显减少，可靠性得到了提升。
%--------------------------------------------------------------------
%9-7a
figure（5）
fortime2=1:
3plot（linspace（0,1/Ts,K*N2）,10*log10（abs（fft（data256（:
time2）,K*N2））.^2/length（data256（:
time2））））;
holdon;
end
holdoff;
xlabel（'frequency'）;
ylabel（'spectraldensityestimate（dB）'）;
xlim（[0（1/Ts/2）]）;
ylim（[-2020]）;
title（'9-7aperiodogram,N=256,rectwindow,nosmooth'）;
%9-7c
figure（7）;
fortime2=1:
3
Sf=abs（fft（data256（:
time2）,K*N2））.^2/length（data256（:
time2））;
forn=1:
（length（Sf）-16）
Sfmat（n,:
）=Sf（n:
n+15,:
）;
end
Sfsmooth=Sfmat*ones（16,1）./16;
plot（linspace（0,1/Ts,n）,10*log10（Sfsmooth））;
holdon;
end
holdoff;
xlabel（'frequency'）;
ylabel（'spectraldensityestimate（dB）'）;
xlim（[0（1/Ts/2）]）;
ylim（[-2020]）;
title（'9-7cperiodogram,N=256,rectwindow,smoothM=2'）;
%--------------------------------------------------------------------
%%9-7b
wn2=2*hanning（N2）;
figure（6）;
fortime2=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;
end
holdoff;
xlabel（'frequency'）;
ylabel（'spectraldensityestimate（dB）'）;
xlim（[0（1/Ts/2）]）;
ylim（[-2020]）;
title（'9-7bperiodogram,N=256,hanningwindow,nosmooth'）;
%9-7d
wn2=2*hanning（N2）;
figure（8）;
fortime2=1:
3
datatapering256（:
time2）=data256（:
time2）.*wn2;Sf=abs（fft（datatapering256（:
time2）,K*N2））.^2/length（datatapering256（:
time2））;
forn=1:
（length（Sf）-16）
Sfmat（n,:
）=Sf（n:
n+15,:
）;
end
Sfsmooth=Sfmat*ones（16,1）./16;
plot（linspace（0,1/Ts,n）,10*log10（Sfsmooth））;
holdon;
end
holdoff;
xlabel（'frequency'）;
ylabel（'spectraldensityestimate（dB）'）;
xlim（[0（1/Ts/2）]）;
ylim（[-2020]）;
title（'9-7dperiodogram,N=256,hanningwindow,smoothM=2'）
从上面四个不同处理的周期图中也可以看到，作升余弦窗处理后的周期图的泄露明显减少。
频率平滑处理提高了可靠性。
在处理数据点数N=64和N=256纵向比较来看，利用的数据越多，谱分析的结果可靠性越高，这符合常理。
%2.sine-waveremoval
clear;clc;clf;
%toloaddata
load（'experiment_data'）;
%data_to_pro=[data（65:
128）data（129:
192）data（193:
256）];t=64;%chosendatalength!
!
!
!
!
data_to_pro=[data（257:
512）data（513:
768）data（769:
1024）];t=256;%chosendatalength!
!
!
!
!
N=8;%paddingnumber
forindex=1:
3;%tocontrolwhichsegmenttouse
freq=（0:
1/（N*t-1）:
1）';
pxx=10*log10（（abs（fft（data_to_pro（:
index）,t*N）））.^2./t）;
figure
（2）
subplot（221）
plot（freq,pxx）;holdon;
axis（[0,0.5,-20,20]）;grid;
title（'PSDoftheoriginaldata'）;
n=（1:
t）';
%-----------------thefirstestimate--------------
[a,i]=max（pxx）;%estimatedfrequsingformula234
pulse=zeros（length（pxx）,1）;pulse（i）=a;istore=ones（length（pxx）,1）;istore（i）=0;%torecordpulseinordertoreinsertafterwards
a=10^（a/10）;%estimatedaisnotindB
fcap=i/（t*N）*2*pi;
sincap=sin（fcap.*n-（t-1）/2*fcap）;%estimatedsinwaveinformula233
coscap=cos（fcap.*n-（t-1）/2*fcap）;%estimatedcoswaveinformula233
pusi=-atan（（sincap'*data_to_pro（:
index））/（coscap'*data_to_pro（:
index）））/pi*180;%estimatephase
acap=2*sqrt（a/t）;%formula235
estimate=acap.*cos（fcap.*n+pusi/180*pi-（t-1）/2*fcap）;%thisisestimatedcosinewaveusea,f,pusi
%--------------thefirstremoval--------------------
removal1_add=data_to_pro（:
index）+estimate;%sincethephaseestimatewillprobablyintroducea'pi'error,whetheraddorsubstracttheestimatesinewaveisunclear
removal1_sub=data_to_pro（:
index）-estimate;
if（removal1_add'*removal1_add>removal1_sub'*removal1_sub）%correctremovalwilldecreasethepowerofvector'data'%
removal1=removal1_sub;%thereforeaftercomparethetwopossibleresultwecandecidewhichonetochose
else
removal1=removal1_add;
end
figure
（2）
subplot（222）
pxx=10*log10（（abs（fft（removal1,t*N）））.^2./t）;
plot（freq,pxx.*istore+pulse）;axis（[0,0.5,-20,20]）;holdon;gridon;
title（'PSDafterthefirstremoval'）
%-----------thesecondestimate----------------
[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）;
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）;
%------------thesecondremoval-------------------------
removal2_add=removal1+estimate;
removal2_sub=removal1-estimate;
if（removal2_add'*removal2_add>removal2_sub'*removal2_sub）
removal2=removal2_sub;
else
removal2=removal2_add;
end
figure
（2）
subplot（223）
pxx=10*log10（（abs（fft（removal2,t*N）））.^2./t）;
plot（freq,pxx.*istore+pulse）;axis（[0,0.5,-20,20]）;holdon;gridon;
title（'PSDaftersecondremoval'）
%-----------------thethirdestimate--------around0.1fs----------
[a,i]=max（pxx（1:
length（pxx）/8））;
pulse（i）=a;istore（i）=0;
a=10^（a/10）;
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_add>removal3_sub'*removal3_sub）
removal3=removal3_sub;
else
removal3=removal3_add;
end
figure
（2）
subplot（224）
pxx=10*log10（（abs（fft（removal3,t*N）））.^2./t）;
plot（freq,pxx.*istore+pulse）;axis（[0,0.5,-20,20]）;holdon;gridon;
title（'9-9aPSDafterthirdremoval'）
end
可以看到用这个方法可以在很大程度上减少泄露。
但是观察N=64的图，可以发现在0.2/T附近的谱线的位置还是很多变，虽然相对之前的方法，我们可以很明显地看到一些分离开来的谱线，但是还是不理想。
N=256的图则效果好很多，0.2/T附近的谱线基本都是一致的，没有很大的变化，从这里也看出了样本的数目对结果还是有很大影响的。
%3.minimum-leakagespectrumestimates
M=16;
f=linspace（0,0.5,1024）;
s=exp（1j*2*pi*（0:
（M-1））'*f）;
fortime1=1:
3
index=N1:
-1:
M;
fortime2=1:
length（index）
X（time2,:
）=data64（index（time2）:
-1:
（index（time2）-M+1）,time1）.';
end
[temp1,temp2]=size（X）;
R=1/temp1*（X'*X）;
Sf（:
time1）=abs（M*1./diag（（s.'*（R\conj（s）））））;
end
figure;
fortime1=1:
3
plot（linspace（0,0.5,1024）,10*log10（Sf（:
time1）））;
holdon;
end
ylim（[-3010]）;
xlim（[00.5]）;
title（'9.10aminimum-leakagespectrumestimateswithN=64M=16'）;
holdoff
比较图9-10a和9-10b，可以看出相同阶数时增大数据量能提高可靠性。
图9-10b~d可以看出滤波器阶数的增加会增大分辨力，但会使可靠性降低，高斯过程的起伏增大，0.1处的谱线峰值降低，这样我们选取滤波器阶数时需要一个折中，使频谱具有相对较好的分辨力和可靠性。
对比谱分析的结果，可以看到minimum-leakagespectrumestimates的优势还是比较明显的，要好于参量方法中的AR建模方法。
%4.ARMethod
functionSf=yw_AR_zty（data,M,f）
%writebyforY-WARMethod
%example£ºdata=data64（:
1）;M=16;f=linspace（0,1,1000）;plot（f,Sf）;
f=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）;
%递推
forN=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）;
forp=1:
N-1
a（N,p）=a（N-1,p）+a（N,N）*a（N-1,N-p）;
end
end
b2（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-11a
figure
（2）;
fortime1=1:
3
M=16;
f=linspace（0,1,1000）;
Sf=yw_AR_zty（data64（:
time1）,M,f）;
plot（f,10*log10（Sf））;
holdon;
end
holdoff;
xlabel（'frequency'）;
ylabel（'spectraldensityestimate（dB）'）;
xlim（[0（1/Ts/2）]）;
ylim（[-2020]）;
title（'9-11aY-WARspectrumestimates,N=64,M=16'）
%9-11b
figure（3）;
fortime2=1:
3
M=16;
f=linspace（0,1,1000）;
Sf=yw_AR_zty（data256（:
time2）,M,f）;
plot（f,10*log10（Sf））;
holdon;
end
holdoff;
xlabel（'frequency'）;
ylab
升级会员 