计算功率谱

clear;
fs=500; %采样率
df=1; %频率分辨率
N=floor(fs/df)+1; %计算的序列点数
t=0:1/fs:(N-1)/fs; %截取信号的时间段
f=0:df:fs; %功率谱估计的频率分辨率和范围
xt=sin(2*pi*50*5)+2*sin(2*pi*130*t)+randn(1,length(t));
%截取时间段上的离散信号样本序列
%
%
%
%利用周期图法进行功率谱估计，但是其得出的功率谱很不平滑，相应的估计协方差比较大
%增加采样点数也不能使周期图变得更加平滑，这是周期图法的缺点，在后面对其改进。
Px=abs(fft(xt)).^2/(N^2); %功率谱估计
Pav_tm=sum(xt.^2)/N; %在时域计算信号功率
Pav_fn=sum(Px); %通过功率谱计算信号功率
figure(1)
subplot(221), %作出功率谱密度图
plot(f,10*log10(Px));xlabel('频率（Hz）');ylabel('功率谱（dB）');
title('周期图法得出的功率谱估计');
%
%
%
%对周期图法进行改进的思想是将信号分段进行估计，然后再将这些估计结果进行平均，
%从而减小估计的协方差，使估计功率谱图变得平滑。
%本程序是将以上501点的信号分为3段，分别作周期图法估计，然后平均。
Px=(abs(fft(xt(1:167))).^2+abs(fft(xt(168:334))).^2+...
abs(fft(xt(335:501))).^2)/3/((N/3)^2);
%分为3段，每段采样点数为N的1/3
Pav_tm=sum(xt.^2)/N; %在时域计算信号功率
Pav_fn=sum(Px); %通过功率谱计算信号功率
subplot(222), %作出功率谱密度图
plot(0:3:fs,10*log10(Px));xlabel('频率（Hz）');ylabel('功率谱（dB）');
title('采用分段估计平均的方法降低估计协方差');
%
%
%
%增加分段数可以进一步降低估计的协方差，然而每段中的数据点太少，就会使估计的
%频率分辨率下降很多。在样本信号序列总点数一定的条件下，可以采用使分段相互重
%叠的方法来增加分段数，从而保持每段中信号点数不变，这样就在保证频率分辨率的
%前提下进一步降低了估计的协方差。本程序中，从分段长度的一半处进行分段重叠，
%这样将501点的信号总共划分为5段，每段167点，相邻段之间重叠83个点。
Px=(abs(fft(xt(1:167))).^2+abs(fft(xt(83:249))).^2+...
abs(fft(xt(168:334))).^2+abs(fft(xt(250:416))).^2+...
abs(fft(xt(335:501))).^2)/5/((N/3)^2);
%看似分为5段，但每段的采样点数都是N的1/3
Pav_tm=sum(xt.^2)/N; %在时域计算信号功率
Pav_fn=sum(Px); %通过功率谱计算信号功率
subplot(223), %作出功率谱密度图
plot(0:3:fs,10*log10(Px));xlabel('频率（Hz）');ylabel('功率谱(dB)');
title('采用重叠分段来增加分段数进一步降低估计协方差');
%
%
%
%采用窗函数对时域信号进

行预处理也可以降低估计的协方差，这也是窗函数的一个应
%用方面。在计算周期图法之前，对数据分段并加非矩形窗（如海明窗等），然后再用
%分段长度一半的混叠率，就能够大大降低估计协方差。这种方法称为Welch平均修正
%周期图法，简称Welch法。本程序利用了汉宁窗。
w=hanning(167)'; %加汉宁窗
w=w*sqrt(167/sum(w.*w)); %使加窗后不对信号功率产生影响
Px=(abs(fft(w.*xt(1:167))).^2+abs(fft(w.*xt(83:249))).^2+...
abs(fft(w.*xt(168:334))).^2+abs(fft(w.*xt(250:416))).^2+...
abs(fft(w.*xt(335:501))).^2)/5/((N/3)^2);
Pav_tm=sum(xt.^2)/N; %在时域计算信号功率
Pav_fn=sum(Px); %通过功率谱计算信号功率
subplot(224); %作出功率谱密度图
plot(0:3:fs,10*log10(Px));xlabel('频率（Hz）');ylabel('功率谱(dB)');
title('采用汉宁窗的welch平均修正周期图法得出的频率谱估计');
%
%
%
%下面的程序是利用MATLAB提供的psd函数和pwelch函数进行功率谱估计
%对信号序列xt采用512点FFT，采样率为500Hz，并使用汉宁窗256点，
%分段混叠点数为128点的welch平均修正周期图法估计频谱
[Px_psd,f_psd]=psd(xt,512,500,hanning(256),128);
%利用psd函数进行频谱估计
[Px_pwelch,f_pwelch]=pwelch(xt,hanning(256),128,512,500);
%利用pwelch函数进行频谱估计
figure(2);
subplot(211), %作出功率谱密度图
plot(f_psd,10*log10(Px_psd/(512/2)));
xlabel('频率（Hz）');ylabel('功率谱(dB)');
title('采用MATLAB的psd指令得出的功率谱估计');
subplot(212),
plot(f_pwelch,10*log10(Px_pwelch/(512/2)));
xlabel('频率（Hz）');ylabel('功率谱(dB)');
title('采用MATLAB的pwelch指令得出的功率谱估计');
%通过功率谱计算信号功率
Pav_psd=sum(Px_psd)/(512/2) %psd函数需要对FFT点数进行归一化处理，即除以FFT点数的一半
Pav_pwelch=sum(Px_pwelch) %pwelch函数已经对FFT点数进行了归一化处理




!!!!!!!!!!!!!!!!!!!
matlab实现经典功率谱估计
1、直接法：
直接法又称周期图法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以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));
window=boxcar(length(xn)); %矩形窗
nfft=1024;
[Pxx,f]=periodogram(xn,window,nfft,Fs); %直接法
plot(f,10*log10(Pxx));

2、间接法：
间接法先由序列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);

3、改进的直接法：
对于直接法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N太小，谱的分辨率又不好，因此需要改进。
3.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_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]);

3.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);