频谱分析

信号的分类

信号分类的方法有很多，这里只给出四种后面将会涉及的分类方式。

按数学关系：确定性信号和随即信号；

按能量特性：能量信号和功率信号；

按取值特征：连续信号和离散信号；

按重复特性，周期性信号和非周期性信号。

随机信号：不能用确定的数学关系来描述的，其值的变动服从统计规律，常用一些数学统计量来描述该类信号的特征。

周期信号：每隔一定时间T，周而复始且无始无终的信号。

以下将重点阐述能量信号和功率信号的概念。我们选择一个具体的物理系统来引出这两种信号的概念。已知阻值为R的电阻上的电压和电流分别为v(t)和i(t)，则通过电阻的瞬时功率为：p(t)=v²(t)/R=i²(t)R。方便起见，假设电阻R为1Ω。故在1Ω电阻上消耗的能量为：

                                                                      

其中f(t)为电压或者电流。

注意积分区间是从，如果积分值存在，则该信号为能量信号；如果该积分不存在，则考虑求该信号的平均功率，如下式：

                                                                  

若该值存在，则该信号即为功率信号。

几点备注：

ü  能量信号的平均功率为0，即E/T（E有限，T无限）为0；

ü  功率信号的能量为无穷，即PT(P有限，T无限)无穷；

ü  信号f(t)可以既非功率信号也非能量信号（如单位斜坡信号f(t)=t）E和P都不存在；

ü  信号f(t)不可能同时既是功率信号也是能量信号（容易理解）；

ü  一般来说，周期信号和随机信号都是功率信号；而非周期的确定信号是能量信号。

ü  除了具有无限能量和功率的信号外，非周期信号要么是能量信号要么是功率信号；

时域和频域

时域和频域是信号通常的两个表示方式。

所谓时域就是自变量是时间，而纵轴是信号的变化，描述了信号在不同时刻取值的函数。

频域的自变量是频率，纵轴是该频率的振幅或相位或功率或能量，描述了一个信号的频谱结构。

如何将信号从时域转换到频域？

一般通过傅里叶变换将信号从时间域转换到频率域。注意，并不是所有的信号都能满足傅里叶变换的条件，通常情况下周期信号满足傅里叶变换的条件。对于测量中经常碰到的随机信号，并不满足傅里叶变换的条件。但是通常会截取其中一段（即样本，用样本的统计值来估计整体信号的频谱特征），而截取的这一段可以进行傅里叶变换，进而可进行功率谱分析。

将信号从时间域转换到频率域后，信号即可表示成很多的谐波之和。

振幅谱：用频域做横坐标，每个谐波的振幅做纵坐标，绘出的图形称为振幅谱；

相位谱：用频域做横坐标，每个谐波的初始相位做纵坐标，绘出的图形称为相位谱。

功率谱和能量谱：

将(1)和(2)式进行傅里叶变换得到信号能量和功率在频域的表示方式，这里不介绍推导过程，只讲述这两个概念的意思。

能量谱密度，也称能量谱，表示信号能量在各个频率的分布情况，信号的总能量等于各个频率分量单独贡献出的能量的连续和。

功率谱密度：也称功率谱，反映了信号功率在各个频率的分布情况，信号的平均功率等于各个频率分量单独贡献出的功率之连续和。

随机信号的特征

随机信号是一类持续时间无限长，具有无限大能量的功率信号，它不满足傅里叶变换条件，而且也不存在解析表达式，因此就不能应用确定信号的频谱计算方法去分析随机信号的频谱。然而，虽然随机信号的频谱不存在，但其相关函数确实是确定的。如果随机信号时平稳的，那么其相关函数的傅里叶变换就是它的功率谱密度函数。

由于在测量中经常用到的是随机信号，这里简单给出描述随机信号的数学指标。这些指标通常包含期望、方差、相关函数、自协方差、互相关函数、互协方差等。

功率谱密度估计

本文以下将详细介绍功率谱估计的方法及其matlab实现。功率谱估计的方法可以分为经典谱估计法和现代谱估计法。

首先，给出一个结论：随机信号自相关函数的傅里叶变换为信号的功率谱密度，互相关函数和互功率谱密度也是一傅里叶变换对。

经典谱估计法概述：

经典谱估计法最常见的有BT(Blackman-Tukey)法和周期图法，它们都是通过傅里叶变换来实现的。BT法是间接法，它先由采样数据估计不同延迟的自相关函数，然后用不同的窗函数对自相关估值开窗加权，再对加权后的自相关估值做傅里叶变换，得到功率谱估计。周期图法是直接法，它是直接将采样数据（也可经窗函数加权）做傅里叶变换，再取其幅度平方而得到功率谱估计。两种方法的结果是一致的。

然而，以上两种方法求出的功率谱分辨率低，故出现了一些改进的方法。大体改进思路是：对数据进行分段，每段求取功率谱密度，继而求平均功率谱密度；数据分段，每段用窗函数进行加权，继而求平均。此外，分割的每段数据还可以有数据重叠。

以上就是经典谱估计方法的发展历程，以下依次用简单的matlab例子进行说明。经典谱估计法的常用matlab函数有periodogram和pwelch函数。

1、 直接法：

直接法又称周期图法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。
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));  %10*log10(Pxx)  %将结果转化为分贝（dB）

为了进一步说明直接法的原理，这里也给出了不用periodogram函数的直接法部分程序：

pxx1=abs(fft(xn(1:256),Nsec).^2)/Nsec;   %根据psd的定义

pxx2=abs(fft(xn(256:512),Nsec).^2)/Nsec;

pxx3=abs(fft(xn(512:768),Nsec).^2)/Nsec;

pxx4=abs(fft(xn(768:1024),Nsec).^2)/Nsec;

pxx=10*log¹⁰((pxx1+pxx2+pxx3+pxx4)/4)  % 转化为分贝

 2、间接法：
间接法先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，便得到x(n)的功率谱估计。
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)分段求周期图再平均。

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]=pwelch(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)，并再周期图计算前直接加进去，加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时，可使各段之间有重叠，这样会使方差减小。

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

分辨率与频谱泄露

在周期图法中，对有限长0-(N-1)的随机数据序列可以看做无限长的随机数据序列经矩形窗开窗截断的结果，我们知道两个时间序列相乘的傅里叶变换等于两个时间序列傅里叶变换的褶积。如果信号的真正功率集中在一个窄的频带内，则该褶积运算将把这个窄带的功率扩展到邻近的频段，这种现象称之为“泄露现象”。弱信号的主瓣很容易被强信号泄露到邻近的副瓣所淹没或畸变。对于短观测数据序列，功率谱轨迹的分辨率是不高的。

结论

在经典谱估计中，总体存在着频率分辨率低、方差性能不好的问题，原因是谱估计时需要对数据进行加窗处理，用有限个数据或其自相关函数来估计无限个数据的功率谱。