信号处理MATLAB实战——功率谱分析
信号处理MATLAB实战——功率谱分析
作者：陈子飞
谨以此文给海洋所于非老师。
转载的同学请保留上面这句话，谢谢。如果还能保留文章来源就更感激不尽了。
1 前言
早先，模拟信道转换成相应的数字信号，工程师们需要知道信号的特征。所以，功率谱分析便由此诞生。
MATLAB中FFT（快速傅里叶变换）是专门用来进行功率谱分析。FFT的基础是DFT（离散傅里叶变换）。在MATLAB中该方法是非常简单的，但是对于许多工程师们来说，在频率与能量计算上仍然存在模糊的认识。
正确地应用功率谱处理信号，首先需要清晰掌握四个概念：频域（frequency domain）、幅值谱（magnitude spectrum）、功率谱（power spectrum）、功率谱密度（power spectrum density）。
2 数学基础
功率谱分析的基础是傅里叶变换：时域到频域的变换
再一个是傅里叶级数，在实数中表示为：
利用欧拉公式，傅里叶级数在复数中的表示：
T为函数的周期，Fn为傅里叶展开系数，
综上可知，时域信号经过傅里叶变换后与时域周期的比值（FFT/T）就是简谐波振幅（很多程序在这里都出现了错误，这也导致谱值的量级出现比较大的差别），即信号的某个简谐波幅值。注意：此时是连续傅里叶变换，下文的离散傅里叶变换也同样存在类似的差异，下文再提。
以上的数学基础就够我们基于MATLAB进行基本操作了。更多的理论与概念可以参考以下链接：
https://baike.baidu.com/item/%E5%82%85%E9%87%8C%E5%8F%B6%E5%8F%98%E6%8D%A2#1
https://blog.csdn.net/liusandian/article/details/51788953
https://www.cnblogs.com/h2zZhou/p/8405717.html
https://www.dsprelated.com/showarticle/1221.php
https://www.dsprelated.com/showarticle/1333.php
https://sites.ualberta.ca/~mlipsett/ENGM541/Power_spectral_density_Matlab.pdf
3 MATLAB应用
接下来，用实例对上面四个概念进行逐一讲解。
3.1频域与幅值谱
clc;clear;
%%
Fs = 1000; % Sampling frequency
T = 1/Fs; % Sampling period
L = 1000; % Length of signal
t = (0:L-1)*T; % Time vector
S = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
X = S + 2*randn(size(t));
subplot(211)
plot(1000*t(1:50),X(1:50))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('t (milliseconds)')
ylabel('X(t)')
%%
Y = fft(S);
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;
subplot(212)
plot(f,P1)
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('|P1(f)|')
（参见：MATLAB help fft example）
图1
MATLAB中实现功率谱分析的基本函数是fft( )，上面程序构建了一个信号，由一个振幅为0.7，频率为50和一个振幅为1,频率为120的两个正弦波组成，并在信号上加入了一个随机噪声，信号在时域上的图像见图1a，可见合成后的信号在时域空间内是无法分辨出是由哪些周期信号组成的，因此，我们对该信号进行频谱分析，见图1b，可见50与120频率的正弦信号在频域中被分解出来。
图1b的y轴含义即幅值，随频域的变化就是幅值谱。这里，提到了两个概念，频域和幅值谱。
这两个概念不难理解，重要的是计算。
提及频域前，先概述一下采样频率fs。假定：时域信号是以最小Nyquist采样，即采样频率是时域信号中最大频率的2倍。根据经验法则，要求信号中研究的频率不得大于采样频率的40%。频域的分辨率δR=fs/nfft。nfft是fft变换时的采样点数。因为fft是基于复平面进行变换，复平面在实数平面投影将出现关于原点对称的幅值谱。所以，幅值谱是关于频域（-0.5 fs: δR : 0.5 fs）对称，所以在取正半轴幅值谱时，幅值需要乘以2。此时的频域为（0：δR ：0.5 fs -δR）。因为0点，所以最后的0.5 fs不包含在频带上。
其次是幅值谱，fft变换是基于离散傅里叶变换DFT进行的，程序中可见，P2 = abs(Y/L)，傅里叶变换的幅值的模与信号点数的比值，准确的说应该是fft变换是采用的点数，即nfft。这是为什么呢？
首先，离散傅里叶变换的定义本身比连续傅里叶变换少了一个dt(采样时间间隔）；
其次，由数学基础部分我们可以知道，连续傅里叶变换时，正弦信号幅值大小与fft变换幅值之间相差一倍周期T，即fft变换的幅值相比实际正弦幅值放大了T倍。所以，通过DFT变换计算幅值谱最终被放大了T/dt=L倍，所以结果需要除以L（做谱的点数）。当然，DFT在选择某些窗口变换时，幅值谱可能还需要考虑恢复系数，这就需要具体问题具体分析了。
综上，我们就清楚了如何计算前两个重要变量——频域和幅值谱。
3.2功率谱与功率谱密度
接下来，继续讲清功率谱与功率谱密度。这段时间，从网上下载的程序可以看出，不少初学者会把这两个概念混为一谈。其实这两者有关联，但也是存在区别的。
功率谱，顾名思义，计算的是功率，或者假若是去除均值后的数据，也可以说是方差。同样，我们采用实例讲解：这是电学领域一个例子，这里，忽略本身物理特性，只看数学运算。
%spectrum_demo.m 1/3/19 Neil Robertson
% Use FFT to find spectrum of sine + noise in units of dBW/bin and dBW/Hz.
fs= 4000; % Hz sample frequency
Ts= 1/fs;
f0= 500; % Hz sine frequency
A= sqrt(2 % V sine amplitude for P= 1 W into 1 ohm.
N= 1024; % number of time samples
nfft= N; % number of frequency samples
n= 0:N-1; % time index
%
x= A*sin(2*pi*f0*n*Ts) + .1*randn(1,N); % 1 W sinewave + noise
w= hanning(N); % window function
window= w.*sqrt(N/sum(w.^2)); % normalize window for P= 1
xw= x.*window'; % apply window to x
%
X= fft(xw,nfft); % DFT
X= X(1:nfft/2); % retain samples from 0 to fs/2
magsq= real(X).^2 + imag(X).^2; % DFT magnitude squared
P_bin= 2/nfft.^2 * magsq; % W/bin power spectrum into 1 ohm
P_Hz= P_bin*nfft/fs; % W/Hz power spectrum
%
PdB_bin= 10*log10(P_bin); % dBW/bin
PdB_Hz= 10*log10(P_Hz); % dBW/Hz
%
k= 0:nfft/2-1; % frequency index
f= k*fs/nfft; % Hz frequency vector
%
plot(f,PdB_bin),grid
axis([0 fs/2 -60 10])
xlabel('Hz'),ylabel('dBW/bin')
title('Spectrum with amplitude units of dBW/bin')
figure
y= max(PdB_Hz);
plot(f,PdB_Hz),grid
axis([0 fs/2 y-60 y+10])
xlabel('Hz'),ylabel('dBW/Hz')
title('Spectrum with amplitude units of dBW/Hz')
程序参见：https://www.dsprelated.com/showarticle/1221.php
开门见山，第一个感叹号强调的P_bin代表功率谱，第二个感叹号强调的P_Hz代表功率谱密度。
可见功率谱与幅值谱的关系为：功率谱=幅值谱模的平方/谱变换点数的平方，在实数空间中，功率谱=2*幅值谱模的平方/谱变换点数的平方，单位为W（watt），反映的是能量。
功率谱密度与功率谱的关系：功率谱密度=功率谱密度*谱变换的点数/采样频率，在实数空间中，功率谱密度=2*功率谱密度*谱变换的点数/采样频率，单位为W/Hz，反映的是单位频率上的能量。
显然，此时我们就清楚了幅值谱、功率谱、功率谱密度三者之间的关系及计算。网上不少程序在这里都出现了错误，比如求功率谱时少除了一个nfft，求功率谱密度时少除了采样频率，这些错误常常发生。
3.3 功率谱密度估计
功率谱密度一般为估计法，因为没有一种方法可以精确计算，现有多种算法估计，包括periodogram, Welch’s method, Yule-Walker Autoregressive Method, Burg Method, etc.。尽管周期图法（基本的fft方法）对大家比较熟悉，但对低信噪比或高方差是，效果也是不理想的（参见https://www.bitweenie.com/listings/power-spectral-density-matlab/）。
MATLAB中常常使用pwelch方法估算功率谱密度。
[pxx,f]= pwelch(x,window,noverlap,nfft,fs);
输入变量:
x input signal vector
window window vector of length nfft
noverlap number of overlapped samples (used for DFT averaging)
nfft number of points in DFT
fs sample rate in Hz
For this article, we’ll assume x is real. Let length(x) = N. If DFT averaging is not desired, you set nfft equal to N and pwelch takes a single DFT of the windowed vector x. If DFT averaging is desired, you set nfft to a fraction of N and pwelch uses windowed segments of x of length nfft as inputs to the DFT. The |X(k)|2 of the multiple DFT’s are then averaged.
输出变量：
pxx power spectral density vector, W/Hz
f vector of frequency values from 0 to fs/2, Hz
下面，我们介绍一下常用的周期图法与pwelch法估算PSD的差异。用一个网络上的例子做分析。需要指出的是，该例子对采样频率、采样点概念上存在一些错误，但不影响算法间的比较结果。实例如下：
周期图法：
clear; clc;
Fs = 1024;
N = Fs*100;
nfft = N;
n = (0:N-1)/N;
Xn = cos(2*pi*100*n)+3*cos(2*pi*200*n)+2*randn(size(n));
figure('color',[1 1 1])
subplot(2,1,1);
plot(Xn);
axis([0 length(Xn) -10 10]);
title('Signal with Noise'); grid on
%%
CXf = fft(Xn, length(Xn));
CXf = abs(CXf);
cpsd = 2*CXf.^2/length(Xn)/Fs;
subplot(2,1,2);
semilogy(cpsd);
axis([0 600 10^(-2) 10^(5)]);
title('Using Periodogram');grid on
Pwelch法
clear; clc;
close all;
Fs = 1024;
N = Fs*100;
nfft = N;
n = (0:N-1)/N;
Xn = cos(2*pi*100*n)+3*cos(2*pi*200*n)+2*randn(size(n));
figure('color',[1 1 1])
subplot(2,1,1);
plot(Xn);
axis([0 length(Xn) -10 10]);
title('Signal with Noise'); grid on
%%
range = 'onesided';
index = round(nfft/2-1);
window_length = 20480;
noverlap = window_length*0.4;
window1 = boxcar(window_length);
[Pxx1,f] = pwelch(Xn,window1,noverlap,nfft,Fs,range);
subplot(212)
semilogy(Pxx1(1:index)); % 10*log10(Pxx)Ϊ·Ö±´´óС£¬´Ë´¦Îª·Ö±´µÄ1/10
axis([0 600 10^(-5) 10^(7)]);
title('Welch-Boxcar');grid on
图3
图4
图3和图4分别是周期图法和pwelch法给出的PSD估计。可见，两种方法均能体现出100hz和200hz的正弦信号。但是，PSD估计值的大小存在差异。上述方法的周期图法仅做了一次周期，增加次数到4次，结果差异明显减小，见图5。
图5
综上，在计算功率谱密度时，推荐使用pwelch方法，目前MATLAB中的pwelch算法相对其他方法较完善，且有多重参数调整方案。
4 结尾
上文均采用谱计算点数等于采样点数计算幅值谱，即nfft=N。当然也可以选用nfft不等于N，此时，fft变换时需要对结果进行平均。定义nfft后，Pwelch函数自动计算DFT平均，这里不再累赘。功率谱分析只是信号处理分析的知识体系的一节，未来，在科研需求的情况下，将继续完善其它章节。
下面，附上一个概念清晰，计算严格的八次DFT变换，求平均的例子，供参考研究：
%spectrum_demo_avg 1/3/19 Neil Robertson
% Perform DFT averaging on signal = sine + noise
% noverlap = 0
% Given signal x(1:8192),
%
% |x(1) x(1025) . . . . x(7169)|
% |x(2) x(1026) . . . . x(7170)|
% xMatrix= | . . . |
% | . . . |
% |x(1024) x(2048) . . . . x(8192)|
%
%
fs= 4000; % Hz sample frequency
Ts= 1/fs;
f0= 500; % Hz sine frequency
A= sqrt(2) % V sine amplitude for P= 1 W into 1 ohm
nfft= 1024; % number of frequency samples
Navg= 8; % number of DFT's to be averaged
N= nfft*Navg; % number of time samples
n= 0:N-1; % time index
x= A*sin(2*pi*f0*n*Ts) + .1*randn(1,N); % 1 W sinewave + noise
w= hanning(nfft); % window function
window= w.*sqrt(nfft/sum(w.^2)); % normalize window for P= 1 W
%
% Create matrix xMatrix with Navg columns,
% each column a segment of x of length nfft
xMatrix= reshape(x,nfft,Navg);
magsq_sum= zeros(nfft/2);
for i= 1:Navg
x_column= xMatrix(:,i);
xw= x_column.*window; % apply window of length nfft
X= fft(xw); % DFT
X= X(1:nfft/2); % retain samples from 0 to fs/2
magsq= real(X).^2 + imag(X).^2; % DFT magnitude squared
magsq_sum= magsq_sum + magsq; % sum of DFT mag squared
end
mag_sq_avg= magsq_sum/Navg; % average of DFT mag squared
P_bin= 2/nfft.^2 *mag_sq_avg; % W/bin power spectrum
P_Hz= P_bin*nfft/fs; % W/Hz power spectrum
PdB_bin= 10*log10(P_bin); % dBW/bin
PdB_Hz= 10*log10(P_Hz); % dBW/Hz
k= 0:nfft/2 -1; % frequency index
f= k*fs/nfft; % Hz frequency vector
%
plot(f,PdB_bin),grid
axis([0 fs/2 -60 10])
xlabel('Hz'),ylabel('dBW/bin')
title('Averaged Spectrum with amplitude units of dBW/bin')
figure
plot(f,PdB_Hz),grid
axis([0 fs/2 -60 10])
xlabel('Hz'),ylabel('dBW/Hz')
title('Averaged Spectrum with amplitude units of dBW/Hz')
程序参见：https://www.dsprelated.com/showarticle/1221.php
作者：czf_po
