功率谱学习及matlab代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
clear;clc;close all;

%% 产生原始信号
f1=5;f2=2;f3=20; f4=100;                           
fs=500;                  			       % 采样率
N=fs*60;                    			   % 采样点数
t = 0 : 1/fs : (N-1)/fs;       			   % 横坐标t坐标精度
% noise_create= randn(1,length(t));        % 产生含有噪声的序列
y1=5*sin(2*pi*f1*t)+10*sin(2*pi*f2*t)+8*sin(2*pi*f3*t)+15*sin(2*pi*f4*t); 
figure(1);subplot(2,1,1);plot(y1,'r');     % 画出原信号的时频图

%% 画出原信号的频率谱
n=0:N-1;
f=n*fs/N;                      			   % 频谱图的频率序列
y=fft(y1,N);                  			   % 对信号进行快速傅里叶变换
mag=abs(y)*2/N;                			   % 求得傅里叶变换后的振幅，与真实幅值相比需要做abs(y)*2/N的运算
subplot(2,1,2); plot(f(1:N/2),mag(1:N/2)); % 绘出Nyquist频率之前随频率变化的振幅

1
2
3
4
5
6
7
8
9
10
%%  直接法(周期图法)求功率谱
% 直接将信号的采样数据x(n)进行Fourier变换求取功率谱密度估计的方法                         
nfft=length(y1);                           % 采样的点数
window=boxcar(length(y1));                 % 使用矩形窗(默认为矩形窗，窗的大小与信号长度一样)
[Pxx,f]=periodogram(y1,window,nfft,fs);    % 直接法
figure(2);subplot(2,2,1);plot(y1);title('信号波形');
subplot(2,2,2);plot(f,10*log10(Pxx));title('周期图法原始信号的功率谱');  %绘制分贝形式的功率谱
% pxx和fft后的平方是一个量纲(w/hz)，10*log10就是用db做单位。

1
2
3
4
5
6
7
8
9
10
11
%%  间接法求功率谱(先估计出自相关函数，再傅里叶变化求功率谱)
nfft=length(y1);
cxn=xcorr(y1,'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)); 
subplot(2,2,3);plot(y1);title('信号波形');
subplot(2,2,4);plot(k,plot_Pxx);title('间接法原始信号的功率谱');     %绘制分贝形式的功率谱

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
%% 改进直接法：Bartlett法（将N点的有限长序列x(n)分段求周期图再平均。）
% 对于直接法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N太小，谱的分辨率又不好，因此需要改进。
nfft=length(y1);
window=boxcar(length(y1));                              %矩形窗
noverlap=0;                                            %数据无重叠
[Pxx,Pxxc]=spectrogram(y1,window,noverlap,nfft,fs);
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(4);
subplot(1,2,1);
plot(k,plot_Pxx);title('直接法信号');
subplot(1,2,2);
plot(k,[plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc]);title('Bartlett信号');

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
%% 加窗平均周期图法 Welch法（选择适当的窗函数，并在周期图计算前直接加进去）
% % Welch法对Bartlett法进行了两方面的修正，一是选择适当的窗函数w(n)w(n)w(n)，并再周期图计算前直接加进去，
% % 加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时，可使各段之间有重叠，这样会使方差减小。
Fs=500;
n=0:1/Fs:1;
nfft=length(y1);
window=boxcar(N);                                         % 矩形窗（时间变量的零次幂窗）
window1=hamming(N);                                       % 海明窗（改进的升余弦窗） 
window2=blackman(N);                                      % blackman窗
noverlap=0;                                               % 数据无重叠
range='onesided';                                         % 频率间隔为[0 Fs/2]，只计算一半的频率
[Pxx,f]=pwelch(y1,window,noverlap,nfft,Fs,range);
[Pxx1,f1]=pwelch(y1,window1,noverlap,nfft,Fs,range);
[Pxx2,f2]=pwelch(y1,window2,noverlap,nfft,Fs,range);
plot_Pxx=10*log10(Pxx);
plot_Pxx1=10*log10(Pxx1);
plot_Pxx2=10*log10(Pxx2);
figure(5);
subplot(1,3,1);
plot(f,plot_Pxx);title('矩形窗');                         % 主瓣比较集中，缺点是旁瓣较高，并有负旁瓣，导致变换中带进了高频干扰和泄漏，甚至出现负谱现象。

subplot(1,3,2);                                           % 海明窗的频谱也是由3个矩形时窗的频谱合成，但其旁瓣衰减速度为20dB/（10oct）
plot(f1,plot_Pxx1);title('汉明窗');                       % (改善频率泄露)加上汉明窗，只有中间的数据体现出来了，两边的数据信息丢失了，所以等会移窗的时候，只会移1/3或1/2窗，这样被前一帧或二帧丢失的数据又重新得到了体现。

subplot(1,3,3);
plot(f2,plot_Pxx2);title('blackman窗');

1
2
3
4
y=fft(y1,N);               
A1=abs(y1);           
A2=A1.^2;