功率谱分析Matlab.docx
《功率谱分析Matlab.docx》由会员分享，可在线阅读，更多相关《功率谱分析Matlab.docx（16页珍藏版）》请在冰豆网上搜索。
功率谱分析Matlab
功率谱分析
由题目内容，设采样频率fs=1000HZ，数据长度为256，模型阶数为14，f1=200,f2=300、250。
（1）用最大熵法进行谱估计
运行程序后，观察图像f1和f2相差较小时，功率谱变化更剧烈；模型的阶数越高，图像中能够获得的信息就越多，但同时计算量也就越大；增加数据长度可以获得更多的信息，提高了谱分析的分辨率，这是因为AR模型的谱估计隐含着对数据和自相关函数的外推，其长度可能会超过给定长度，分辨率不受信源信号的限制。
（2）分别用Levinson递推法和Burg法进行功率谱分析
Levinson递推法
运行程序后，观察图像，f1和f2相差较小时，功率谱变化更剧烈；模型的阶数越高，图像中能够获得的信息就越多，但同时计算量也就越大；增加数据长度可以获得更多的信息，提高了谱分析的分辨率，但本题中信号为正弦信号加白噪声，故图像观察不明显。
Burg法
运行程序后，观察图像，f1和f2相差较小时，功率谱变化更剧烈；模型的阶数越高，图像中能够获得的信息就越多，但同时计算量也就越大；增加数据长度可以获得更多的信息，提高了谱分析的分辨率。
（3）改变信号的相位、频率、信噪比，上述谱分析结果有何变化
如果正弦信号的频率过大，超过fs/2，会产生频率混叠现象，输入f1=600HZ，会在400HZ处产生一个波峰；降低信噪比会导致谱分辨率下降；信号起始相位的变动可导致谱线的偏移和分裂（我的图像观察不到）。
最大熵法估计
N=1024;
Nfft=256;
Fs=1000;
n=0:
N-1;
t=n/Fs;
x1=sin（2*pi*200*t）;
x2=sin（2*pi*300*t）;%0.3
xn=x1+awgn（x1,10）+x2+awgn（x2,10）;
[Pxx1,f]=pmem（xn,14,Nfft,Fs）;
subplot（4,1,1）
plot（f,10*log10（Pxx1））;
xlabel（'Frequency（Hz）'）;ylabel（'PowerSpectrum（dB）'）;
title（'MEMf2/fs=0.3,Nfft=256,Oder=14'）;
grid
N=1024;
Nfft=256;
Fs=1000;
n=0:
N-1;
t=n/Fs;
x1=sin（2*pi*200*t）;
x2=sin（2*pi*250*t）;%0.25
xn=x1+awgn（x1,10）+x2+awgn（x2,10）;
[Pxx1,f]=pmem（xn,14,Nfft,Fs）;
subplot（4,1,2）
plot（f,10*log10（Pxx1））;
xlabel（'Frequency（Hz）'）;ylabel（'PowerSpectrum（dB）'）;
title（'MEMf2/fs=0.25,Nfft=256,Oder=14'）;
grid
N=1024;
Nfft=512;%修改数据长度512
Fs=1000;
n=0:
N-1;
t=n/Fs;
x1=sin（2*pi*200*t）;
x2=sin（2*pi*300*t）;%0.3
xn=x1+awgn（x1,10）+x2+awgn（x2,10）;
[Pxx1,f]=pmem（xn,14,Nfft,Fs）;
subplot（4,1,3）
plot（f,10*log10（Pxx1））;
xlabel（'Frequency（Hz）'）;ylabel（'PowerSpectrum（dB）'）;
title（'MEMf2/fs=0.3,Nfft=512,Oder=14'）;
grid
N=1024;
Nfft=256;
Fs=1000;
n=0:
N-1;
t=n/Fs;
x1=sin（2*pi*200*t）;
x2=sin（2*pi*300*t）;%0.3
xn=x1+awgn（x1,10）+x2+awgn（x2,10）;
[Pxx1,f]=pmem（xn,24,Nfft,Fs）;%修改阶数为24
subplot（4,1,4）
plot（f,10*log10（Pxx1））;
xlabel（'Frequency（Hz）'）;ylabel（'PowerSpectrum（dB）'）;
title（'MEMf2/fs=0.3,Nfft=256,Oder=24'）;
Grid
Burg法估计
N=1024;
Nfft=256;
Fs=1000;
n=0:
N-1;
t=n/Fs;
x1=sin（2*pi*200*t）;
x2=sin（2*pi*300*t）;%0.3
xn=x1+awgn（x1,10）+x2+awgn（x2,10）;
[Pxx1,f]=pburg（xn,14,Nfft,Fs）;
subplot（4,1,1）
plot（f,10*log10（Pxx1））;
xlabel（'Frequency（Hz）'）;ylabel（'PowerSpectrum（dB）'）;
title（'Burgf2/fs=300,Nfft=256,Oder=14'）;
grid
N=1024;
Nfft=256;
Fs=1000;
n=0:
N-1;
t=n/Fs;
x1=sin（2*pi*200*t）;
x2=sin（2*pi*250*t）;%0.25
xn=x1+awgn（x1,10）+x2+awgn（x2,10）;
[Pxx1,f]=pburg（xn,14,Nfft,Fs）;
subplot（4,1,2）
plot（f,10*log10（Pxx1））;
xlabel（'Frequency（Hz）'）;ylabel（'PowerSpectrum（dB）'）;
title（'Burgf2/fs=250,Nfft=256,Oder=14'）;
grid
N=1024;
Nfft=512;%修改数据长度512
Fs=1000;
n=0:
N-1;
t=n/Fs;
x1=sin（2*pi*200*t）;
x2=sin（2*pi*300*t）;%0.3
xn=x1+awgn（x1,10）+x2+awgn（x2,10）;
[Pxx1,f]=pburg（xn,14,Nfft,Fs）;
subplot（4,1,3）
plot（f,10*log10（Pxx1））;
xlabel（'Frequency（Hz）'）;ylabel（'PowerSpectrum（dB）'）;
title（'Burgf2/fs=300,Nfft=512,Oder=14'）;
grid
N=1024;
Nfft=256;
Fs=1000;
n=0:
N-1;
t=n/Fs;
x1=sin（2*pi*200*t）;
x2=sin（2*pi*300*t）;%0.3
xn=x1+awgn（x1,10）+x2+awgn（x2,10）;
[Pxx1,f]=pburg（xn,24,Nfft,Fs）;%修改阶数为24
subplot（4,1,4）
plot（f,10*log10（Pxx1））;
xlabel（'Frequency（Hz）'）;ylabel（'PowerSpectrum（dB）'）;
title（'Burgf2/fs=300,Nfft=256,Oder=24'）;
grid
Levinson递推法
N=1024;
Nfft=256;
Fs=1000;
n=0:
N-1;
t=n/Fs;
x1=sin（2*pi*200*t）;
x2=sin（2*pi*300*t）;%0.3
xn=x1+awgn（x1,10）+x2+awgn（x2,10）;
[Pxx1,f]=pyulear（xn,14,Nfft,Fs）;
%[Pxx1,f]=Levinson（xn,14,Nfft,Fs）;
subplot（4,1,1）
plot（f,10*log10（Pxx1））;
xlabel（'Frequency（Hz）'）;ylabel（'PowerSpectrum（dB）'）;
title（'LevinsonNfft=256,f2/fs=0.3,Oder=14'）;
grid
N=1024;
Nfft=256;
Fs=1000;
n=0:
N-1;
t=n/Fs;
x1=sin（2*pi*200*t）;
x2=sin（2*pi*250*t）;%0.25
xn=x1+awgn（x1,10）+x2+awgn（x2,10）;
[Pxx1,f]=pyulear（xn,14,Nfft,Fs）;
subplot（4,1,2）
plot（f,10*log10（Pxx1））;
xlabel（'Frequency（Hz）'）;ylabel（'PowerSpectrum（dB）'）;
title（'LevinsonNfft=256,f2/fs=0.25,Oder=14'）;
grid
N=1024;
Nfft=512;%修改数据长度512
Fs=1000;
n=0:
N-1;
t=n/Fs;
x1=sin（2*pi*200*t）;
x2=sin（2*pi*300*t）;%0.3
xn=x1+awgn（x1,10）+x2+awgn（x2,10）;
[Pxx1,f]=pyulear（xn,14,Nfft,Fs）;
subplot（4,1,3）
plot（f,10*log10（Pxx1））;
xlabel（'Frequency（Hz）'）;ylabel（'PowerSpectrum（dB）'）;
title（'LevinsonNfft=512,f2/fs=0.3,Oder=14'）;
grid
N=1024;
Nfft=256;
Fs=1000;
n=0:
N-1;
t=n/Fs;
x1=sin（2*pi*200*t）;
x2=sin（2*pi*300*t）;%0.3
xn=x1+awgn（x1,10）+x2+awgn（x2,10）;
[Pxx1,f]=pyulear（xn,24,Nfft,Fs）;%修改阶数为24
subplot（4,1,4）
plot（f,10*log10（Pxx1））;
xlabel（'Frequency（Hz）'）;ylabel（'PowerSpectrum（dB）'）;
title（'LevinsonNfft=256,f2/fs=0.3,Oder=24'）;
grid
最大熵法改变信号的相位、频率、信噪比
N=1024;
Nfft=256;
Fs=1000;
n=0:
N-1;
t=n/Fs;
x1=sin（2*pi*200*t）;
x2=sin（2*pi*300*t）;%0.3
xn=x1+awgn（x1,10）+x2+awgn（x2,10）;
[Pxx1,f]=pmem（xn,14,Nfft,Fs）;
subplot（4,1,1）
plot（f,10*log10（Pxx1））;
xlabel（'Frequency（Hz）'）;ylabel（'PowerSpectrum（dB）'）;
title（'MEMf2/fs=0.3,Nfft=256,Oder=14'）;
grid
N=1024;
Nfft=256;
Fs=1000;
n=0:
N-1;
t=n/Fs;
x1=sin（2*pi*200*t+pi/6）;%相位加了pi/6
x2=sin（2*pi*300*t）;%0.3
xn=x1+awgn（x1,10）+x2+awgn（x2,10）;
[Pxx1,f]=pmem（xn,14,Nfft,Fs）;
subplot（4,1,2）
plot（f,10*log10（Pxx1））;
xlabel（'Frequency（Hz）'）;ylabel（'PowerSpectrum（dB）'）;
title（'MEMf2/fs=0.3,Nfft=256,Oder=14，相位加pi/6'）;
grid
N=1024;
Nfft=256;
Fs=1000;
n=0:
N-1;
t=n/Fs;
x1=sin（2*pi*200*t）;
x2=sin（2*pi*300*t）;%0.3
xn=x1+awgn（x1,5）+x2+awgn（x2,5）;%性噪比改为5
[Pxx1,f]=pmem（xn,14,Nfft,Fs）;
subplot（4,1,3）
plot（f,10*log10（Pxx1））;
xlabel（'Frequency（Hz）'）;ylabel（'PowerSpectrum（dB）'）;
title（'MEMf2/fs=0.3,Nfft=256,Oder=14，性噪比=5'）;
grid
N=1024;
Nfft=256;
Fs=1000;
n=0:
N-1;
t=n/Fs;
x1=sin（2*pi*300*t）;
x2=sin（2*pi*400*t）;%0.3
xn=x1+awgn（x1,10）+x2+awgn（x2,10）;
[Pxx1,f]=pmem（xn,14,Nfft,Fs）;
subplot（4,1,4）
plot（f,10*log10（Pxx1））;
xlabel（'Frequency（Hz）'）;ylabel（'PowerSpectrum（dB）'）;
title（'MEMf1/fs=0.3,f2/fs=0.4,Nfft=256,Oder=14'）;
grid
Burg改变信号的相位、频率、信噪比
N=1024;
Nfft=256;
Fs=1000;
n=0:
N-1;
t=n/Fs;
x1=sin（2*pi*200*t）;
x2=sin（2*pi*300*t）;%0.3
xn=x1+awgn（x1,10）+x2+awgn（x2,10）;
[Pxx1,f]=pburg（xn,14,Nfft,Fs）;
subplot（4,1,1）
plot（f,10*log10（Pxx1））;
xlabel（'Frequency（Hz）'）;ylabel（'PowerSpectrum（dB）'）;
title（'Burgf2/fs=300,Nfft=256,Oder=14'）;
grid
N=1024;
Nfft=256;
Fs=1000;
n=0:
N-1;
t=n/Fs;
x1=sin（2*pi*200*t+pi/6）;%相位加了pi/6
x2=sin（2*pi*300*t）;%0.3
xn=x1+awgn（x1,10）+x2+awgn（x2,10）;
[Pxx1,f]=pburg（xn,14,Nfft,Fs）;
subplot（4,1,2）
plot（f,10*log10（Pxx1））;
xlabel（'Frequency（Hz）'）;ylabel（'PowerSpectrum（dB）'）;
title（'Burgf2/fs=300,Nfft=256,Oder=14，相位加pi/6'）;
grid
N=1024;
Nfft=256;
Fs=1000;
n=0:
N-1;
t=n/Fs;
x1=sin（2*pi*300*t）;
x2=sin（2*pi*400*t）;%0.3
xn=x1+awgn（x1,10）+x2+awgn（x2,10）;
[Pxx1,f]=pburg（xn,14,Nfft,Fs）;
subplot（4,1,3）
plot（f,10*log10（Pxx1））;
xlabel（'Frequency（Hz）'）;ylabel（'PowerSpectrum（dB）'）;
title（'Burgf1/fs=300,f2/fs=400,Nfft=256,Oder=14'）;
grid
N=1024;
Nfft=256;
Fs=1000;
n=0:
N-1;
t=n/Fs;
x1=sin（2*pi*200*t）;
x2=sin（2*pi*300*t）;%0.3
xn=x1+awgn（x1,5）+x2+awgn（x2,5）;%性噪比改为5
[Pxx1,f]=pburg（xn,14,Nfft,Fs）;
subplot（4,1,4）
plot（f,10*log10（Pxx1））;
xlabel（'Frequency（Hz）'）;ylabel（'PowerSpectrum（dB）'）;
title（'Burgf2/fs=300,Nfft=256,Oder=14,性噪比=5'）;
grid
Levinson法改变信号的相位、频率、信噪比
N=1024;
Nfft=256;
Fs=1000;
n=0:
N-1;
t=n/Fs;
x1=sin（2*pi*200*t）;
x2=sin（2*pi*300*t）;%0.3
xn=x1+awgn（x1,10）+x2+awgn（x2,10）;
[Pxx1,f]=pyulear（xn,14,Nfft,Fs）;
%[Pxx1,f]=Levinson（xn,14,Nfft,Fs）;
subplot（4,1,1）
plot（f,10*log10（Pxx1））;
xlabel（'Frequency（Hz）'）;ylabel（'PowerSpectrum（dB）'）;
title（'LevinsonNfft=256,f2/fs=0.3,Oder=14'）;
grid
N=1024;
Nfft=256;
Fs=1000;
n=0:
N-1;
t=n/Fs;
x1=sin（2*pi*200*t+pi/6）;%相位加了pi/6
x2=sin（2*pi*300*t）;%0.3
xn=x1+awgn（x1,10）+x2+awgn（x2,10）;
[Pxx1,f]=pyulear（xn,14,Nfft,Fs）;
%[Pxx1,f]=Levinson（xn,14,Nfft,Fs）;
subplot（4,1,2）
plot（f,10*log10（Pxx1））;
xlabel（'Frequency（Hz）'）;ylabel（'PowerSpectrum（dB）'）;
title（'LevinsonNfft=256,f2/fs=0.3,Oder=14，相位加pi/6'）;
grid
N=1024;
Nfft=256;
Fs=1000;
n=0:
N-1;
t=n/Fs;
x1=sin（2*pi*300*t）;
x2=sin（2*pi*400*t）;%0.3
xn=x1+awgn（x1,10）+x2+awgn（x2,10）;
[Pxx1,f]=pyulear（xn,14,Nfft,Fs）;
%[Pxx1,f]=Levinson（xn,14,Nfft,Fs）;
subplot（4,1,3）
plot（f,10*log10（Pxx1））;
xlabel（'Frequency（Hz）'）;ylabel（'PowerSpectrum（dB）'）;
title（'LevinsonNfft=256,f1/fs=0.3,f2/fs=0.4,Oder=14'）;
grid
N=1024;
Nfft=256;
Fs=1000;
n=0:
N-1;
t=n/Fs;
x1=sin（2*pi*200*t）;
x2=sin（2*pi*300*t）;%0.3
xn=x1+awgn（x1,5）+x2+awgn（x2,5）;%性噪比位5
[Pxx1,f]=pyulear（xn,14,Nfft,Fs）;
%[Pxx1,f]=Levinson（xn,14,Nfft,Fs）;
subplot（4,1,4）
plot（f,10*log10（Pxx1））;
xlabel（'Frequency（Hz）'）;ylabel（'PowerSpectrum（dB）'）;
title（'LevinsonNfft=256,f2/fs=0.3,Oder=14,性噪比=5'）;
grid
升级会员 