随机信号利用经典谱估计法估计信号的功率谱作业综述：给出一段信号“asd.wav”,利用经典谱估计法的原理，通过不同的谱估计方法，求出信号的功率谱密度函数。

采用MATLAB语言，利用MATLAB语言强大的数据处理和数据可视化能力，通过GUI的对话框模板，使操作更为简便！在一个GUI界面中，同时呈现出不同方法产生出的功率谱。

这里给出了几种不同的方法：BT法，周期图法，平均法以及Welch法。

把几种不同方法所得到的功率谱都呈现在一个界面中，便于对几种不同方法得到的功率谱作对比。

一．题目要求给出一段信号及采样率，利用经典谱估计法估计出信号的功率谱。

二．基本原理及方法经典谱估计的方法，实质上依赖于传统的傅里叶变换法。

它是将数据工作区外的未知数据假设为零，相当于数据加窗，主要方法有BT法，周期图法，平均法以及Welch法。

1. BT法(Blackman-Tukey)●理论基础：(1)随机序列的维纳-辛钦定理由于随机序列{X(n)}的自相关函数Rx(m)=E[X(n)X(n+m)]定义在离散点m上，设取样间隔为，则可将随机序列的自相关函数用连续时间函数表示为等式两边取傅里叶变换，则随机序列的功率谱密度(2)谱估计BT法是先估计自相关函数Rx(m)(|m|=0,1,2…,N-1),然后再经过离散傅里叶变换求的功率谱密度的估值。

即其中可有式得到。

2. 周期图法●理论基础：周期图法是根据各态历经随机过程功率谱的定义来进行谱估计的。

在前面我们已知，各态历经的连续随机过程的功率谱密度满足式中 是连续随机过程第i 个样本的截取函数 的频谱。

对应在随机序列中则有由于随机序列中观测数据 仅在 的点上存在，则 的N 点离散傅里叶变换为：因此有随机信号的观测数据 的功率谱估计值（称“周期图”）如下：由于上式中的离散傅里叶变换可以用快速傅里叶变换计算，因此就可以估计出功率谱。

3.平均法：理论基础：平均法可视为周期图法的改进。

周期图经过平均后会使它的方差减少，达到一致估计的目的，有一个定理：如果 ， ， ， 是不相关的随机变量，且都有个均值 及其方差 ，则可以证明它们的算术平均的均值为 ，方差为。

由定理可见：具有 个独立同分布随机变量平均的方差，是单个随机变量方差的 ，当 时，方差，可以达到一致估计的目的。

因此，将 个独立的估计量经过算术平均后得到的估计量的方差也是原估计量方差的 。

平均图法即是将数据 ， ， 分段求周期图法后再平均。

例如，给定N=1000个数据样本（平均法适用于数据量大的场合），则可以将它分成10个长度为100的小段，分别计算每一段的周期图()()21001100,100(1)1,1,2,```,10100l j l n l G w X e l ω-=-==∑然后将这10个周期图加以平均得谱估计值：()()10100100,1110l l G w G w ==∑由于这10小段的周期图取决于同一个过程，因而其均值相同。

若这10个小段的周期图是统计独立的，则这10个小段平均之后的方差却是单段方差的 。

D[ ]= D[ ]即：平均法将 ， ， 的N 个数据分成L 段（N=ML ），若各数据段相互独立，则平方后估计量的方差是原来不分段估计量方差的 。

所以当 时，估计量的方差趋于0，达到一致估计的目的。

但是，随着分段数L 的增加，M 点数减少，分辨率减少，使估计变成有偏估计。

相反，若L 减少，M 增加，虽偏差减少，但方差增大。

所以，在应用中，必须兼顾分辨率和方差的要求来适当选择M 和L 的值。

4.Welch 法：理论基础：Welch 法又称修正周期法，其步骤为：○1先将N 个数据分成L 段，每数据段M 个数据，N=ML ； ○2选择适当的窗函数w （n ），并用该w （n ）依次对每段数据做相应的加权，然后确定每段的周期图()21,(1)1(),1,2,```,Ml j n M l n n M l G w x w n e l L MUω--=-==∑U 为归一化因子1201(),M n U w n M-==∑对每段周期图进行平均得到功率谱估计：()(),11LM M l l G w G w L ==∑当数据量一定时，若分段数L 增加，M 点数减少，则分辨率下降；若L 减少，虽M 增加，但方差增大。

解决这一矛盾的方法是，让数据段间适当“重叠”。

三．算法设计与实现1程序流程图（1） BT 法的流程图（2）周期图法的流程图（3）平均法的流程图（4）Welch法的流程图2.主要模块的设计:(1)产生原始信号Fs=600;nfft=1024;n=0:1/Fs:1;y=cos(2*pi*30*n)+3*cos(2*pi*100*n)+randn(size(n));axes(handles.axes1);plot(n,y);title(‘原始信号波形‘);原始信号波形如下：(2)BT法global y Fs nfft n;cxn=xcorr(y,'unbiased');CXk=fft(cxn,nfft);Pxx=abs(CXk);index=0:round(nfft/2-1);k=index*Fs/nfft;plot_Pxx=10*log(Pxx(index+1))/(log(10)); axes(handles.axes2);plot(k,plot_Pxx);xlabel('频率（hz）');ylabel('功率谱密度（Db/Hz)');title('BT法');BT法所得功率谱波形如下：(3)周期图法global y Fs nfft n;XF=fft(y,nfft);Pxx=abs(XF).^2/length(n);index=0:round(nfft/2-1);f=index*Fs/nfft;axes(handles.axes4);plot(f,10*log(Pxx(index+1)));xlabel('频率（hz）');ylabel('功率谱密度（Db/Hz)');title('周期图法');周期图法所得功率谱波形如下：(4)平均法global y Fs nfft n;window=hamming(nfft);noverlap=0;p=0.9;[Pxx,Pxxc]=psd(y,nfft,Fs,window,noverlap,p); index=0:round(nfft/2-1);k=index*Fs/nfft;Qxx=10*log10(Pxx(index+1));axes(handles.axes6);plot(k,Qxx );xlabel('频率（hz）');ylabel('功率谱密度（Db/Hz)');title('平均法');平均法所得功率谱波形如下：(5) Welch法global y Fs nfft n;window1=hamming(100);noverlap=20;range='half';[Pxx1,f]=pwelch(y,window1,noverlap,nfft,Fs,range);Qxx1=10*log10(Pxx1);axes(handles.axes5);plot(f,Qxx1);xlabel('频率（hz）');ylabel('功率谱密度（Db/Hz)');title('Welch法');Welch法所得功率谱波形如下：四. 软件使用说明综述从给出一段信号y=cos(2*pi*30*n)+3*cos(2*pi*100*n),利用经典谱估计法的原理，通过不同的谱估计方法，求出信号的功率谱密度函数，并利用GUI界面呈现出不同谱估计方法所得的结果。

BT法BT法是先估计自相关函数Rx(m)(|m|=0,1,2…,N-1),然后再经过离散傅里叶变换求的功率谱密度的估值。

由式=估算出，再对作FFT变换，得到。

周期图法周期图法是根据各态历经随机过程功率谱的定义来进行谱估计的。

获取x(n)后，对x(n)作FFT得到X(W)，再由得出功率谱。

平均法平均法可视为周期图法的改进。

周期图经过平均后会使它的方差减少，达到一致估计的目的。

获取x(n)后，将x(n)分为10段，对每段用 … 计算出周期图，对以上10个周期图加以平均得出功率谱。

Welch 法Welch 法又称修正周期法。

获取x(n)后，先将N 个数据分成L 段，选择汉明窗w(n )，并用该w(n )依次对每段数据做相应的加权，确定每段的周期图()21,(1)1(),1,2,```,Ml j n M l n n M l G w x w n e l L MU ω--=-==∑由()(),11LM M l l G w G w L ==∑得出每段谱估计。

GUI 界面如下：GUI 界面五. 结果分析由图可以看出，在频率30hz 和100hz 处功率谱有两个峰值，说明信号中有30hz 和100hz 的周期成分。

通过BT法能观察到两个峰值，但是所呈现的波形不能准确表达出信号的功率谱变化情况。

通过周期图法求出的功率谱密度在很大范围内波动，而且容易证明，即使增加信号取样点数N，实验效果依然没有明显改进。

用有限长样本序列的DFT来表示随机序列的功率谱只是一种估计或近似，不可避免存在误差。

为了减少误差，使功率谱估计更加平滑，可采用平均法。

平均法采用了分段的功率谱估计，较之于周期图法，平均法的估计曲线较为平滑。

在程序中加入窗函数，使得谱分辨率不会由于分段而下降。

Welch法就是用改进的平均法来求取随机信号的功率谱密度估计。

Welch法采用信号重叠分段，加窗函数和FFT算法等计算一个信号序列的功率谱。

六. 任务分工王龙元（学号：3832008009）完成了BT法以及周期图法的程序编写以及GUI界面的设计和PPT的编写;郭敏（学号：3222008054）完成了平均法以及Welch法的程序编写以及实验报告的统筹,编写.七. 参考文献[1] 常建平, 李林海. 随机信号分析. 科学出版社.2010年，第5章第3节：184－189；[2] 黄文梅, 熊桂林, 杨勇.信号分析与处理.国防科技大学出版社.2000年,第6章第3节:222—229八．附录程序源代码：function varargout = ming(varargin)gui_Singleton = 1;gui_State = struct('gui_Name', mfilename, ...'gui_Singleton', gui_Singleton, ...'gui_OpeningFcn', @ming_OpeningFcn, ...'gui_OutputFcn', @ming_OutputFcn, ...'gui_LayoutFcn', [] , ...'gui_Callback', []);if nargin && ischar(varargin{1})gui_State.gui_Callback = str2func(varargin{1});endif nargout[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});elsegui_mainfcn(gui_State, varargin{:});endfunction ming_OpeningFcn(hObject, eventdata, handles, varargin) guidata(hObject, handles);function varargout = ming_OutputFcn(hObject, eventdata, handles) varargout{1} = handles.output;function figure1_CreateFcn(hObject, eventdata, handles)function pushbutton1_Callback(hObject, eventdata, handles)global y Fs nfft n;Fs=600;nfft=1024;n=0:1/Fs:1;y=cos(2*pi*30*n)+3*cos(2*pi*100*n)+randn(size(n)); axes(handles.axes1);plot(n,y);title(‘原始信号波形‘);function pushbutton3_Callback(hObject, eventdata, handles)global y Fs nfft n;window=hamming(nfft);noverlap=0;p=0.9;[Pxx,Pxxc]=psd(y,nfft,Fs,window,noverlap,p);index=0:round(nfft/2-1);k=index*Fs/nfft;Qxx=10*log10(Pxx(index+1));axes(handles.axes6);plot(k,Qxx );xlabel('频率（hz）');ylabel('功率谱密度（Db/Hz)');title('平均法');function pushbutton5_Callback(hObject, eventdata, handles)global y Fs nfft n;cxn=xcorr(y,'unbiased');CXk=fft(cxn,nfft);Pxx=abs(CXk);index=0:round(nfft/2-1);k=index*Fs/nfft;plot_Pxx=10*log(Pxx(index+1))/(log(10));axes(handles.axes2);plot(k,plot_Pxx);xlabel('频率（hz）');ylabel('功率谱密度（Db/Hz)');title('BT法');function pushbutton4_Callback(hObject, eventdata, handles) global y Fs nfft n;window1=hamming(100);noverlap=20;range='half';[Pxx1,f]=pwelch(y,window1,noverlap,nfft,Fs,range);Qxx1=10*log10(Pxx1);axes(handles.axes5);plot(f,Qxx1);xlabel('频率（hz）');ylabel('功率谱密度（Db/Hz)');title('Welch法');function pushbutton6_Callback(hObject, eventdata, handles) global y Fs nfft n;XF=fft(y,nfft);Pxx=abs(XF).^2/length(n);index=0:round(nfft/2-1);f=index*Fs/nfft;axes(handles.axes4);plot(f,10*log(Pxx(index+1)));xlabel('频率（hz）');ylabel('功率谱密度（Db/Hz)');title('周期图法');欢迎您的下载，资料仅供参考！致力为企业和个人提供合同协议，策划案计划书，学习资料等等打造全网一站式需求。