matlab的功率谱计算
功率谱估计在现代信号处理中是⼀个很重要的课题，涉及的问题很多。在这⾥，结合matlab，我做⼀个粗略介绍。功率谱估计可以分为经典谱估计⽅法与现代谱估计⽅法。经典谱估计中最简单的就是周期图法，⼜分为直接法与间接法。直接法先取N点数据的傅⾥叶变换（即频谱），然后取频谱与其共轭的乘积，就得到功率谱的估计；间接法先计算N点样本数据的⾃相关函数，然后取⾃相关函数的傅⾥叶变换，即得到功率谱的估计.都可以编程实现，很简单。在matlab中，周期图法可以⽤函数periodogram实现。
但是周期图法估计出的功率谱不够精细，分辨率⽐较低。因此需要对周期图法进⾏修正，可以将信号序列x(n)分为n个不相重叠的⼩段，分别⽤周期图法进⾏谱估计，然后将这n段数据估计的结果的平均值作为整段数据功率谱估计的结果。还可以将信号序列x(n)重叠分段，分别计算功率谱，再计算平均值作为整段数据的功率谱估计。
这2种称为分段平均周期图法，⼀般后者⽐前者效果好。加窗平均周期图法是对分段平均周期图法的改进，即在数据分段后，对每段数据加⼀个⾮矩形窗进⾏预处理，然后在按分段平均周期图法估计功率谱。相对于分段平均周期图法，加窗平均周期图法可以减⼩频率泄漏，增加频峰的宽度。welch法就是利⽤改进的平均周期图法估计估计随机信号的功率谱，它采⽤信号分段重叠，加窗，FFT等技术来计算功
率谱。与周期图法⽐较，welch法可以改善估计谱曲线的光滑性，⼤⼤提⾼谱估计的分辨率。matlab中，welch法⽤函数psd实现。调⽤格式如下：[Pxx,F] = PSD(X,NFFT,Fs,WINDOW,NOVERLAP)
X：输⼊样本数据
NFFT：FFT点数
Fs：采样率
WINDOW：窗类型
NOVERLAP，重叠长度
现代谱估计主要针对经典谱估计分辨率低和⽅差性不好提出的，可以极⼤的提⾼估计的分辨率和平滑性。可以分为参数模型谱估计和⾮参数模型谱估计。参数模型谱估计有AR模型，MA模型，ARMA模型等；⾮参数模型谱估计有最⼩⽅差法和MUSIC法等。由于涉及的问题太多，这⾥不再详述，可以参考有关资料。
matlab中，现代谱估计的很多⽅法都可以实现。music⽅法⽤pmusic命令实现；pburg函数利⽤burg法实现功率谱估计；pyulear函数利⽤yule-walker算法实现功率谱估计等等。
另外，sptool⼯具箱也具有功率谱估计的功能。窗⼝化的操作界⾯很⽅便，⽽且有多种⽅法可以选择。
t=linspace(0,1,1000);
signal=4*sin(2*pi*50*t)+5*sin(2*pi*200*t);
noise=randn(size(t));
symbol=signal+noise;
Y=fft(symbol,128);
f=1000*(0:64)/128;
P1=Y.*conj(Y)/128; %直接法
subplot(1,3,1);
plot(f,10*log10(P1(1:65)));
xlabel('frequency')
ylabel('power')
title('直接法')
直接法：
直接法⼜称周期图法，它是把随机序列x(n)的N个观测数据视为⼀能量有限的序列，直接计算x(n)的离散傅⽴叶变换，得X(k)，然后再取其幅值的平⽅，并除以N，作为序列x(n)真实功率谱的估计。
Matlab代码⽰例：
clear;
Fs=1000; %采样频率
n=0:1/Fs:1;
%产⽣含有噪声的序列
matlab求傅里叶变换xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
window=boxcar(length(xn)); %矩形窗
nfft=1024;
[Pxx,f]=periodogram(xn,window,nfft,Fs); %直接法
plot(f,10*log10(Pxx));
%间接法
cxn=xcorr(symbol,'unbiased'); %计算序列的⾃相关函数
P2=fft(cxn,128);
subplot(1,3,2);
plot(f,10*log10(P2(1:65)));
xlabel('frequency')
ylabel('power')
title('间接法')
%加窗法
window2=blackman(100); %blackman窗
noverlap=20; %数据⽆重叠
range='onesided'; %频率间隔为[0 1000/2]，只计算⼀半的频率
[P3(1:65),f]=pwelch(symbol,window2,noverlap,128,1000,range);
plot_P3=10*log10(P3(1:65));
subplot(1,3,3);
plot(f,plot_P3(1:65));
xlabel('frequency')
ylabel('power')
title('加窗法')
间接法：
间接法先由序列x(n)估计出⾃相关函数R(n)，然后对R(n)进⾏傅⽴叶变换，便得到x(n)的功率谱估计。
Matlab代码⽰例：
clear;
Fs=1000; %采样频率
n=0:1/Fs:1;
%产⽣含有噪声的序列
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
nfft=1024;
cxn=xcorr(xn,'unbiased'); %计算序列的⾃相关函数
CXk=fft(cxn,nfft);
Pxx=abs(CXk);
index=0:round(nfft/2-1);
k=index*Fs/nfft;
plot_Pxx=10*log10(Pxx(index+1));
plot(k,plot_Pxx);
改进的直接法：
对于直接法的功率谱估计，当数据长度N太⼤时，谱曲线起伏加剧，若N太⼩，谱的分辨率⼜不好，因此需要改进。
1. Bartlett法
Bartlett平均周期图的⽅法是将N点的有限长序列x(n)分段求周期图再平均。
Matlab代码⽰例：
clear；
Fs=1000;
n=0:1/Fs:1;
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
nfft=1024;
window=boxcar(length(n)); %矩形窗
noverlap=0; %数据⽆重叠
p=0.9; %置信概率
[Pxx,Pxxc]=psd(xn,nfft,Fs,window,noverlap,p);
index=0:round(nfft/2-1);
k=index*Fs/nfft;
plot_Pxx=10*log10(Pxx(index+1));
plot_Pxx=10*log10(Pxx(index+1));
plot_Pxxc=10*log10(Pxxc(index+1));
figure(1)
plot(k,plot_Pxx);
pause;
figure(2)
plot(k,[plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc]);
2. Welch法
Welch法对Bartlett法进⾏了两⽅⾯的修正，⼀是选择适当的窗函数w(n)，并再周期图计算前直接加进去，加窗的优点是⽆论什么样的窗函数均可使谱估计⾮负。⼆是在分段时，可使各段之间有重叠，这样会使⽅差减⼩。
Matlab代码⽰例：
clear;
Fs=1000;
n=0:1/Fs:1;
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
nfft=1024;
window=boxcar(100); %矩形窗
window1=hamming(100); %海明窗
window2=blackman(100); %blackman窗
noverlap=20; %数据⽆重叠
range='half'; %频率间隔为[0 Fs/2]，只计算⼀半的频率
[Pxx,f]=pwelch(xn,window,noverlap,nfft,Fs,range);
[Pxx1,f]=pwelch(xn,window1,noverlap,nfft,Fs,range);
[Pxx2,f]=pwelch(xn,window2,noverlap,nfft,Fs,range);
plot_Pxx=10*log10(Pxx);
plot_Pxx1=10*log10(Pxx1);
plot_Pxx2=10*log10(Pxx2);
figure(1)
plot(f,plot_Pxx);
pause;
figure(2)
plot(f,plot_Pxx1);
pause;
figure(3)
plot(f,plot_Pxx2);
版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。
发表评论