《功率谱图应用》由会员分享，可在线阅读，更多相关《功率谱图应用（7页珍藏版）》请在金锄头文库上搜索。
1、1. . 基本方法基本方法 周期图法是直接将信号的采样数据 x(n)进行 Fourier 变换求取功率谱密度估计的方法。 假定有限长随机信号序列为 x(n)。它的 Fourier 变换和功率谱密度估计存在下面的关系：式中，N 为随机信号序列 x(n)的长度。在离散的频率点 f=kf,有：其中，FFTx(n)为对序列 x(n)的 Fourier 变换，由于 FFTx(n)的周期为 N，求得的功率谱估计以 N 为周期，因此这种方法称为周期图法。下面用例子说明如何采用这种方法进行功率谱用有限长样本序列的 Fourier 变换来表示随机序列的功率谱，只是一种估计或近似，不可避免存在误差。为了减少误差，
2、使功率谱估计更加平滑，可采用分段平均周期图法（Bartlett 法）、加窗平均周期图法（Welch 法）等方法加以改进。2. 分段平均周期图法（分段平均周期图法（Bartlett 法）法） 将信号序列 x(n)，n=0,1,N-1，分成互不重叠的 P 个小段，每小段由 m 个采样值，则 P*m=N。对每个小段信号序列进行功率谱估计，然后再取平均作为整个序列 x(n)的功率谱 估计。 平均周期图法还可以对信号 x(n)进行重叠分段，如按 2:1 重叠分段，即前一段信号和后一 段信号有一半是重叠的。对每一小段信号序列进行功率谱估计，然后再取平均值作为整个 序列 x(n)的功率谱估计。这两种方法都称
3、为平均周期图法，一般后者比前者好。程序运行 结果为图 9-5，上图采用不重叠分段法的功率谱估计，下图为 2:1 重叠分段的功率谱估计， 可见后者估计曲线较为平滑。与上例比较，平均周期图法功率谱估计具有明显效果（涨落 曲线靠近 0dB）。 3.3.加窗平均周期图法加窗平均周期图法 加窗平均周期图法是对分段平均周期图法的改进。在信号序列 x(n)分段后，用非矩形窗口 对每一小段信号序列进行预处理，再采用前述分段平均周期图法进行整个信号序列 x(n)的 功率谱估计。由窗函数的基本知识（第 7 章）可知，采用合适的非矩形窗口对信号进行处 理可减小“频谱泄露”，同时可增加频峰的宽度，从而提高频谱分辨率。
4、 其中上图采用无重叠数据分段的加窗平均周期图法进行功率谱估计，而下图采用重叠数据 分段的加窗平均周期图法进行功率谱估计，显然后者是更佳的，信号谱峰加宽，而噪声谱 均在 0dB 附近，更为平坦（注意采用无重叠数据分段噪声的最大的下降分贝数大于 5dB， 而重叠数据分段周期图法噪声的最大下降分贝数小于 5dB）。 4. Welch 法估计及其法估计及其 MATLAB 函数函数 Welch 功率谱密度就是用改进的平均周期图法来求取随机信号的功率谱密度估计的。Welch 法采用信号重叠分段、加窗函数和 FFT 算法等计算一个信号序列的自功率谱估计(PSD 如上 例中的下半部分的求法)和两个信号序列的互
5、功率谱估计（CSD）。 MATLAB 信号处理工具箱函数提供了专门的函数 PSD 和 CSD 自动实现 Welch 法估计，而不需 要自己编程。 （1） 函数 psd 利用 Welch 法估计一个信号自功率谱密度，函数调用格式为：Pxx,f=psd(x,Nfft,Fs,window,Noverlap,dflag) 式中，x 为信号序列；Nfft 为采用的 FFT 长度。这一值决定了功率谱估计速度，当 Nfft 采 用 2 的幂时，程序采用快速算法；Fs 为采样频率；Window 定义窗函数和 x 分段序列的长度。窗函数长度必须小于或等于 Nfft，否则会给出错误信息；Noverlap 为分段序
6、列重叠的采样 点数（长度），它应小于 Nfft;dflag 为去除信号趋势分量的选择项：linear,去除线 性趋势分量，mean去除均值分量，none不做去除趋势处理。Pxx 为信号 x 的自功 率谱密度估计。f 为返回的频率向量，它和 Pxx 对应，并且有相同长度。 在 psd 函数调用格式中，缺省值为：Nfft=min(256,length(x), Fs=2Hz, window=hanning(Nfft),noverlap=0. 若 x 是实序列，函数 psd 仅计算频率为正的 功率 注意程序前半部分中频率向量 f 的创建方法。它与函数 psd 的输出 Pxx 长度的关系如 下：若 x
7、为实序列，当 Nfft 为奇数时，f=(0:(Nfft+1)/2-1)/Nfft;当 Nfft 为偶数时， f=(0:Nfft/2)/Nfft。 函数还有一种缺省返回值的调用格式，用于直接绘制信号序列 x 的功率谱估计曲线。 函数还可以计算带有置信区间的功率谱估计，调用格式为：Pxx,Pxxc,f=psd(x,Nfft,Fs,window,Noverlap,p) 式中，p 为置信区间，0=p=1。由此可知，滤波器输入白噪声序列的输出信号的功率谱或自相关可以确定滤波器的频 率特性。 （2）函数 csd 利用 welch 法估计两个信号的互功率谱密度，函数调用格式为：Pxy,f=csd(x,y,N
8、fft,Fs,window,Noverlap,dflag)Pxy,Pxyc,f=csd(x,y,Nfft,Fs,window,Noverlap,p) 这里，x,y 为两个信号序列；Pxy 为 x,y 的互功率谱估计；其他参数的意义同自功率谱函数 psd。 可以看到，两个白噪声信号的互功率谱（上图）杂乱无章，看不出周期成分，大部分功率 谱在-5dB 以下。然而白噪声与带有噪声的周期信号的功率谱在其周期（频率为 1000Hz）处 有一峰值，清楚地表明了周期信号的周期或频率。因此，利用未知信号与白噪声信号的互 功率谱也可以检测未知信号中所含有的频率成分。5 5 多多 窗窗 口口 法法 多窗口法(Mu
9、ltitaper method,简称 MTM 法)利用多个正交窗口（Tapers）获得各自独立的 近似功率谱估计，然后综合这些估计得到一个序列的功率谱估计。相对于普通的周期图法， 这种功率谱估计具有更大的自由度，并在估计精度和估计波动方面均有较好的效果。普通 的功率谱估计只利用单一窗口，因此在序列始端和末端均会丢失相关信息，而且无法找回。 而 MTM 法估计增加窗口使得丢失的信息尽量减少。 MTM 法简单地采用一个参数：时间带宽积(Time-bandwidth product)NW，这个参数用 以定义计算功率谱所用窗的数目，为 2*NW-1。NW 越大，功率谱计算次数越多，时间域分辨 率越高，
10、而频率域分辨率降低，使得功率谱估计的波动减小。随着 NW 增大，每次估计中谱 泄漏增多，总功率谱估计的偏差增大。对于每一个数据组，通常有一个最优的 NW 使得在估 计偏差和估计波动两方面求得折中，这需要在程序中反复调试来获得。 MATLAB 信号处理工具箱中函数 PMTM 就是采用 MTM 法估计功率谱密度。函数调用格式 为：Pxx,f=pmtm(x,nw,Nfft,Fs) 式中，x 为信号序列；nw 为时间带宽积，缺省值为 4。通常可取 2，5/2，3，7/2；Nfft 为 FFT 长度；Fs 为采样频率。上面的函数还可以通过无返回值而绘出置信区间，如 pmtm(x,nw,Nfft,Fs,o
11、ption,p) 绘制带置信区间的功率谱密度估计曲线，0=p=1。 6 6 最最 大大 熵熵 法（法（MaxmumMaxmum entropyentropy methodmethod， MEMMEM 法）法） 如上所述，周期图法功率谱估计需要对信号序列“截断”或加窗处理，其结果是使估计的 功率谱密度为信号序列真实谱和窗谱的卷积，导致误差的产生。 最大熵功率谱估计的目的是最大限度地保留截断后丢失的“窗口”以外信号的信息，使估 计谱的熵最大。主要方法是以已知的自相关序列 rxx(0),rxx(1),rxx(p)为基础，外推自 相关序列 rxx(p+1),rxx(p+2),，保证信息熵最大。 最大熵
12、功率谱估计法假定随机过程是平稳高斯过程，可以证明，随机信号的最大熵谱与 AR 自回归（全极点滤波器）模型谱是等价的。 MATLAB 信号处理工具箱提供最大熵功率谱估计函数 pmem，其调用格式为：Pxx,f,a=pmem(x,p,Nfft,Fs,xcorr) 式中，x 为输入信号序列或输入相关矩阵；p 为全极点滤波器阶次；a 为全极点滤波器模型 系数向量；xcorr是把 x 认为是相关矩阵。比较最大熵功率谱估计(MEM)和改进的平均周期图功率谱估计，可见，MEM 法估计的功 率谱曲线较光滑。在这一方法中，MEM 法选定全极点滤波器的阶数取得越大，能够获得 的窗口外的信息越多，但计算量也越大，需
13、要根据情况折中考虑。7 7 多信号分类法多信号分类法 MATLAB 信号处理工具箱还提供另一种功率谱估计函数 pmusic。该函数执行多信号分类法 （multiple signal classification， Music 法)。将数据自相关矩阵看成由信号自相关 矩阵和噪声自相关矩阵两部分组成，即数据自相关矩阵 R 包含有两个子空间信息：信号子 空间和噪声子空间。这样，矩阵特征值向量（Eigen vector）也可分为两个子空间：信号 子空间和噪声子空间。为了求得功率谱估计，函数 pmusic 计算信号子空间和噪声子空间的 特征值向量函数，使得在周期信号频率处函数值最大，功率谱估计出现峰值，
14、而在其他频 率处函数值最小。其调用格式为：Pxx,f,a=pmusic(x,p,thresh,Nfft,Fs,window,Noverlap) 式中,x 为输入信号的向量或矩阵；p 为信号子空间维数；thresh 为阈值，其他参数的意义 与函数 psd 相同。功率谱密度相关方法的 MATLAB 实现% %分段平均周期图法（Bartlett 法）%运用信号不重叠分段估计功率谱 Nsec=256;n=0:sigLength-1;t=n/Fs; %数据点数，分段间隔，时间序列pxx1=abs(fft(y(1:256),Nsec).2)/Nsec; %第一段功率谱pxx2=abs(fft(y(257:
15、512),Nsec).2)/Nsec; %第二段功率谱pxx3=abs(fft(y(515:768),Nsec).2)/Nsec; %第三段功率谱pxx4=abs(fft(y(769:1024),Nsec).2)/Nsec; %第四段功率谱Pxx=10*log10(pxx1+pxx2+pxx3+pxx4)/4); %平均得到整个序列功率谱f=(0:length(Pxx)-1)*Fs/length(Pxx); %给出功率谱对应的频率%plot(f(1:Nsec/2),Pxx(1:Nsec/2); %绘制功率谱曲线figure,plot(f(1:Nsec),Pxx(1:Nsec); %绘制功率谱曲
16、线xlabel(频率/Hz);ylabel(功率谱 /dB);title(平均周期图(无重叠) N=4*256);grid on% %运用信号重叠分段估计功率谱 pxx1=abs(fft(y(1:256),Nsec).2)/Nsec; %第一段功率谱pxx2=abs(fft(y(129:384),Nsec).2)/Nsec; %第二段功率谱pxx3=abs(fft(y(257:512),Nsec).2)/Nsec; %第三段功率谱pxx4=abs(fft(y(385:640),Nsec).2)/Nsec; %第四段功率谱pxx5=abs(fft(y(513:768),Nsec).2)/Nsec; %第五段功率谱pxx6=abs(fft(y(641:896),Nsec).2)/Nsec; %第六段功率谱pxx7=abs(fft(y(769:1024)
