功率谱分析

一、直接法
clear;clc;close all;  %清除变量；清屏；关闭当前图形窗口
Fs=1000;
t=0:1/Fs:1;
nfft=2048;  %改变nfft的值可对比不同采样值时的谱估计效果
%****************生成信号、噪声**************%
x1=cos(2*pi*40*t)+3*cos(2*pi*45*t);%信号
x2=randn(size(t));                  %噪声
x3=x1+x2;         %信号+噪声
[Pxx,f]=periodogram(x3,window,nfft,Fs); %直接法
plot(f,10*log10(Pxx));title('直接法 nfft=2048');;
set(gca,'xlim',[1 120]); ;ylabel('Am/dB');
xlabel('Frequency/Hz');
二、间接法
Fs=1000;% 采样频率
n=0:1/Fs:1;% 产生含有噪声的序列
x1=cos(2*pi*40*n)+3*cos(2*pi*45*n);%信号
x2=randn(size(n));                  %噪声
x3=x1+x2;                           %信号+噪声
nfft=1024;
cxn=xcorr(x3);% 计算序列的自相关函数
CXk=fft(cxn);
Pxx=abs(CXk);
index=0:round(nfft/2-1);
f=index*Fs/nfft;
plot_Pxx=10*log10(Pxx(index+1));
figure (1)
plot(f,plot_Pxx);
title('间接法 nfft=1024');ylabel('Am/dB');
set(gca,'xlim',[1 120]);
xlabel('Frequency/Hz');
三、Bartlett法
clear;clc;close all;  %清除变量；清屏；关闭当前图形窗口
Fs=1000;
t=0:1/Fs:1;
nfft=1024;
%****************生成信号、噪声**************%
x1=cos(2*pi*40*t)+3*cos(2*pi*45*t);%信号
x2=randn(size(t));                  %噪声
x3=x1+x2;                           %信号+噪声
window=hamming(512);                    %海明窗
noverlap=0; %数据无重叠
p=0.9; %置信概率
[Pxx,Pxxc]=psd(x3,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);title('Bartlett法海明窗');;
set(gca,'xlim',[1 120]); ;ylabel('Am/dB');
xlabel('Frequency/Hz');
四、Welch法
clear;clc;close all;  %清除变量；清屏；关闭当前图形窗口
Fs=1000;
t=0:1/Fs:1;
nfft=1024;
%****************生成信号、噪声**************%
x1=cos(2*pi*40*t)+3*cos(2*pi*45*t);%信号
x2=randn(size(t));                  %噪声
x3=x1+x2;                           %信号+噪声
window=hamming(512);                    %海明窗
noverlap=128;
range='onesided';                   %频率间隔为[0 Fs/2]，只计算一半的频
[Pxx1,f]=pwelch(x3,window,noverlap,nfft,Fs,range);
plot_Pxx1=10*log10(Pxx1);
figure(1);
plot(f,plot_Pxx1);title('Welch法海明窗');ylabel('Am/dB');
set(gca,'xlim',[1 120]); xlabel('Frequency/Hz');

对所给的实测信号进行谱估计，本文采用了周期图法和Welch法。其程序如下所示：

一、周期图法
clear;clc;close all;%清除变量；清屏；关闭当前图形窗口
Fs=2000; %采样频率
load Chanel8Xia2Data.mat;
x1=Chanel8Xia2Gray;
x2=Chanel8Xia2Sti;
Mlag=length(x1);
%****************求功率谱密度**************%
[Pxx1,f]=periodogram(x1,window,length(x1),Fs);%直接法
[Pxx2,f]=periodogram(x2,window,length(x1),Fs); 
plot_Pxx1=10*log10(Pxx1);
plot_Pxx2=10*log10(Pxx2);
%****************显示功率谱密度曲线**************%
figure(1);
plot(f,plot_Pxx1,'b');
axis([0,100,-50,40]);grid on;
title('直接法 ');;
set(gca,'xlim',[1 100]);
ylabel('Am/dB');
hold on;
plot(f,plot_Pxx2,'r');
axis([0,100,-50,40]);grid on;
title('直接法 ');
xlabel('Frequency/Hz');ylabel('Am/dB'); 
二、Welch法
clear;clc;close all;  %清除变量；清屏；关闭当前图形窗口
Fs=2000;
load Chanel8Xia2Data.mat;
x1=Chanel8Xia2Gray;
x2=Chanel8Xia2Sti;
window1=hamming(1024);                    %海明窗
noverlap=256;                             %数据无重叠
range='onesided';                        %频率间隔为[0 Fs/2]，只计算一半的频、
%****************求功率谱密度**************%
[Pxx1,f]=pwelch(x1,window1,noverlap,length(x1),Fs,range);
[Pxx2,f]=pwelch(x2,window1,noverlap,length(x1),Fs,range);
plot_Pxx1=10*log10(Pxx1);
plot_Pxx2=10*log10(Pxx2);
%****************显示功率谱密度曲线**************%
figure(1);
plot(f,plot_Pxx1,'b');
axis([0,100,-30,30]);grid on;
title('Welch法海明窗');ylabel('Am/dB');
set(gca,'xlim',[1 100]); 
hold on;
plot(f,plot_Pxx2,'r');
title('Welch法海明窗');ylabel('Am/dB');
xlabel('Frequency/Hz');