功率谱密度机器实现
1. 基本方法
周期图法是直接将信号的采样数据x(n)进行Fourier变换求取功率谱密度估计的方法。假定有限长随机信号序列为x(n)。它的Fourier变换和功率谱密度估计存在下面的关系：
式中，N为随机信号序列x(n)的长度。在离散的频率点f=kΔf,有：
其中，FFT[x(n)]为对序列x(n)的Fourier变换，由于FFT[x(n)]的周期为N，求得的功率谱估计以N为周期，因此这种方法称为周期图法。下面用例子说明如何采用这种方法进行功率谱
用有限长样本序列的Fourier变换来表示随机序列的功率谱，只是一种估计或近似，不可避免存在误差。为了减少误差，使功率谱估计更加平滑，可采用分段平均周期图法（Bartlett法）、加窗平均周期图法（Welch 法）等方法加以改进。
2. 分段平均周期图法（Bartlett法）
将信号序列x(n)，n=0,1,…,N-1，分成互不重叠的P个小段，每小段由m个采样值，则P*m=N。对每个小段信号序列进行功率谱估计，然后再取平均作为整个序列x(n)的功率谱估计。
平均周期图法还可以对信号x(n)进行重叠分段，如按2:1重叠分段，即前一段信号和后一段信号有一半是重叠的。对每一小段信号序列进行功率谱估计，然后再取平均值作为整个序列x(n)的功率谱估计。这两种方法都称为平均周期图法，一般后者比前者好。程序运行结果为图9-5，上图采用不重叠分段法的功率谱估计，下图为2:1重叠分段的功率谱估计，可见后者估计曲线较为平滑。与上例比较，平均周期图法功率谱估计具有明显效果（涨落曲线靠近0dB）。
3.加窗平均周期图法
加窗平均周期图法是对分段平均周期图法的改进。在信号序列x(n)分段后，用非矩形窗口对每一小段信号序列进行预处理，再采用前述分段平均周期图法进行整个信号序列x(n)的功率谱估计。由窗函数的基本知识（第7章）可知，采用合适的非矩形窗口对信号进行处理可减小“频谱泄露”，同时可增加频峰的宽度，从而提高频谱分辨率。
其中上图采用无重叠数据分段的加窗平均周期图法进行功率谱估计，而下图采用重叠数据分段的加窗平均周期图法进行功率谱估计，显然后者是更佳的，信号谱峰加宽，而噪声谱均在0dB附近，更为平坦（注意采用无重叠数据分段噪声的最大的下降分贝数大于5dB，而重叠数据分段周期图法噪声的最大下降分贝数小于5dB）。
4. Welch法估计及其MATLAB函数
Welch功率谱密度就是用改进的平均周期图法来求取随机信号的功率谱密度估计的。Welch 法采用信号重叠分段、加窗函数和FFT算法等计算一个信号序列的自功率谱估计(PSD如上例中的下半部分的求法)和两个信号序列的互功率谱估计（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为分段序列重叠的采样
点数（长度），它应小于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为实序列，当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,Nfft,Fs,window,Noverlap,’dflag’)
[Pxy,Pxyc[,f]]=csd(x,y,Nfft,Fs,window,Noverlap,p)
这里，x,y为两个信号序列；Pxy为x,y的互功率谱估计；其他参数的意义同自功率谱函数psd。
可以看到，两个白噪声信号的互功率谱（上图）杂乱无章，看不出周期成分，大部分功率谱在-5dB以下。然而白噪声与带有噪声的周期信号的功率谱在其周期（频率为1000Hz）处有一峰值，清楚地表明了周期信号的周期或频率。因此，利用未知信号与白噪声信号的互功率谱也可以检测未知信号中所含有的频率成分。
5 多窗口法
多窗口法(Multitaper method,简称MTM法)利用多个正交窗口（Tapers）获得各自独立的近似功率谱估计，然后综合这些估计得到一个序列的功率谱估计。相对于普通的周期图法，这种功率谱估计具有更大的自由度，并在估计精度和估计波动方面均有较好的效果。普通的功率谱估计只利用单一窗口，因此在序列始端和末端均会丢失相关信息，而且无法找回。而MTM法估计增加窗口使得丢失的信息尽量减少。
MTM法简单地采用一个参数：时间带宽积(Time-bandwidth product)NW，这个参数用以定义计算功率谱所用窗的数目，为2*NW-1。NW越大，功率谱计算次数越多，时间域分辨率越高，而频率域分辨率降低，使得功率谱估计的波动减小。随着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,’option’,p)绘制带置信区间的功率谱密度估计曲线，0<=p<=1。
6 最大熵法（Maxmum entropy method，MEM法）
如上所述，周期图法功率谱估计需要对信号序列“截断”或加窗处理，其结果是使估计的功率谱密度为信号序列真实谱和窗谱的卷积，导致误差的产生。
最大熵功率谱估计的目的是最大限度地保留截断后丢失的“窗口”以外信号的信息，使估计谱的熵最大。主要方法是以已知的自相关序列r xx(0),r xx(1),…,r xx(p)为基础，外推自相关序列r xx(p+1),r xx(p+2),…，保证信息熵最大。
最大熵功率谱估计法假定随机过程是平稳高斯过程，可以证明，随机信号的最大熵谱与AR 自回归（全极点滤波器）模型谱是等价的。
MATLAB信号处理工具箱提供最大熵功率谱估计函数pmem，其调用格式为：
[Pxx,f,a]=pmem(x,p,Nfft,Fs,’xcorr’)
式中，x为输入信号序列或输入相关矩阵；p为全极点滤波器阶次；a为全极点滤波器模型系数向量；’xcorr’是把x认为是相关矩阵。
比较最大熵功率谱估计(MEM)和改进的平均周期图功率谱估计，可见，MEM法估计的功率谱曲线较光滑。在这一方法中，MEM法选定全极点滤波器的阶数取得越大，能够获得的窗口外的信息越多，但计算量也越大，需要根据情况折中考虑。
7 多信号分类法
MATLAB信号处理工具箱还提供另一种功率谱估计函数pmusic。该函数执行多信号分类法（multiple signal classification，Music法)。将数据自相关矩阵看成由信号自相关矩阵和噪声自相关矩阵两部分组成，即数据自相关矩阵R包含有两个子空间信息：信号子空间和噪声子空间。这样，矩阵特征值向量（Eigen vector）也可分为两个子空间：信号子空间和噪声子空间。为了求得功率谱估计，函数pmusic计算信号子空间和噪声子空间的特征值向量函数，使得在周期信号频率处函数值最大，功率谱估计出现峰值，而在其他频率处函数值最小。其调用格式为：
[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: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)); %绘制功率谱曲线
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),Nsec).^2)/Nsec; %第七段功率谱
Pxx=10*log10((pxx1+pxx2+pxx3+pxx4+pxx5+pxx6+pxx7)/7); %功率谱平均并转化为dB f=(0:length(Pxx)-1)*Fs/length(Pxx); %频率序列
figure,plot(f(1:Nsec/2),Pxx(1:Nsec/2)); %绘制功率谱曲线
xlabel('频率/Hz'); ylabel('功率谱/dB');
title('平均周期图(重叠一半) N=1024');
grid on
%%
Nsec=256;n=0:sigLength-1;t=n/Fs; %数据长度，分段数据长度、时间序列
w=hanning(256); %采用的窗口数据
%采用不重叠加窗方法的功率谱估计
pxx1=abs(fft(w.*y(1:256),Nsec).^2)/norm(w)^2; %第一段加窗振幅谱平方
pxx2=abs(fft(w.*y(257:512),Nsec).^2)/norm(w)^2; %第二段加窗振幅谱平方
pxx3=abs(fft(w.*y(513:768),Nsec).^2)/norm(w)^2; %第三段加窗振幅谱平方
pxx4=abs(fft(w.*y(769:1024),Nsec).^2)/norm(w)^2; %第四段加窗振幅谱平方
Pxx=10*log10((pxx1+pxx2+pxx3+pxx4)/4); %求得平均功率谱，转换为dB
f=(0:length(Pxx)-1)*Fs/length(Pxx); %求得频率序列
figure
subplot(2,1,1),plot(f(1:Nsec/2),Pxx(1:Nsec/2)); %绘制功率谱曲线
xlabel('频率/Hz');ylabel('功率谱/dB');
title('加窗平均周期图(无重叠) N=4*256');
grid on
%采用重叠加窗方法的功率谱估计
pxx1=abs(fft(w.*y(1:256),Nsec).^2)/norm(w)^2; %第一段加窗振幅谱平方
pxx2=abs(fft(w.*y(129:384),Nsec).^2)/norm(w)^2; %第二段加窗振幅谱平方
pxx3=abs(fft(w.*y(257:512),Nsec).^2)/norm(w)^2; %第三段加窗振幅谱平方
pxx4=abs(fft(w.*y(385:640),Nsec).^2)/norm(w)^2; %第四段加窗振幅谱平方
pxx5=abs(fft(w.*y(513:768),Nsec).^2)/norm(w)^2; %第五段加窗振幅谱平方
pxx6=abs(fft(w.*y(641:896),Nsec).^2)/norm(w)^2; %第六段加窗振幅谱平方
pxx7=abs(fft(w.*y(769:1024),Nsec).^2)/norm(w)^2; %第七段加窗振幅谱平方
Pxx=10*log10((pxx1+pxx2+pxx3+pxx4+pxx5+pxx6+pxx7)/7);%平均功率谱转换为dB
f=(0:length(Pxx)-1)*Fs/length(Pxx); %频率序列
subplot(2,1,2),plot(f(1:Nsec/2),Pxx(1:Nsec/2)); %绘制功率谱曲线
xlabel('频率/Hz');ylabel('功率谱/dB');
title('加窗平均周期图(重叠一半)N=1024');
grid on
%%
%4分段平均周期图法（hanning窗）
Nsec=256;
n=0:sigLength-1;
t=n/Fs;
w=hanning(256);
Pxx1=abs(fft(w.*y(1:256),Nsec).^2)/Nsec;
Pxx2=abs(fft(w.*y(257:512),Nsec).^2)/Nsec;
Pxx3=abs(fft(w.*y(513:768),Nsec).^2)/Nsec;
Pxx4=abs(fft(w.*y(769:1024),Nsec).^2)/Nsec;
Pxx=10*log10((Pxx1+Pxx2+Pxx3+Pxx4)/4);
f=(0:length(Pxx)-1)*Fs/length(Pxx);
figure
subplot(2,1,1)
plot(f,Pxx);
xlabel('频率/Hz');ylabel('功率谱/dB');
title('Averaged Modified Periodogram (none overlap) N=4*256'); grid
%4分段（2：1重叠）平均周期图法（hanning窗）
Nsec=256;
n=0:sigLength-1;
t=n/Fs;
w=hanning(256);
Pxx1=abs(fft(w.*y(1:256),Nsec).^2)/Nsec;
Pxx2=abs(fft(w.*y(129:384),Nsec).^2)/Nsec;
Pxx3=abs(fft(w.*y(257:512),Nsec).^2)/Nsec;
Pxx4=abs(fft(w.*y(385:640),Nsec).^2)/Nsec;
Pxx5=abs(fft(w.*y(513:768),Nsec).^2)/Nsec;
Pxx6=abs(fft(w.*y(641:896),Nsec).^2)/Nsec;
Pxx7=abs(fft(w.*y(769:1024),Nsec).^2)/Nsec;
Pxx=10*log10((Pxx1+Pxx2+Pxx3+Pxx4+Pxx5+Pxx6+Pxx7)/7);
f=(0:length(Pxx)-1)*Fs/length(Pxx);
subplot(2,1,2);plot(f,Pxx);
xlabel('频率/Hz');ylabel('Power Spectrum (dB)');
title('Averaged Modified Periodogram (half overlap) N=1024'); grid
%%
%PSD_WELCH方法
%采样频率
Nfft=256;n=0:sigLength-1;t=n/Fs; %数据长度、时间序列
window=hanning(256); %选用的窗口
noverlap=128; %分段序列重叠的采样点数（长度）
dflag='none'; %不做趋势处理
[Pxx,Pxxc,f]=psd(y,Nfft,Fs,window,noverlap,; %功率谱估计,并以的置信度给出置信区间，无返回值是绘制出置信区间
figure
plot(f,10*log10(Pxx)); %绘制功率谱
xlabel('频率/Hz');ylabel('功率谱/dB');
title('PSD—Welch方法'); grid on
%%
%最大熵法（MEM法）
Nfft=256;n=0:sigLength-1;t=n/Fs; %数据长度、分段长度和时间序列
window=hanning(256); %采用窗口
[Pxx1,f]=pmem(x,20,Nfft,Fs); %采用最大熵法，采用滤波器阶数14，估计功率谱figure,subplot(2,1,1),plot(f,10*log10(Pxx1)); %绘制功率谱
xlabel('频率/Hz');ylabel('功率谱/dB');
title('最大熵法Order=20原始信号功率谱');grid on
[Pxx1,f]=pmem(y0,20,Nfft,Fs); %采用最大熵法，采用滤波器阶数14，估计功率谱subplot(2,1,2),plot(f,10*log10(Pxx1)); %绘制功率谱
xlabel('频率/Hz');ylabel('功率谱/dB');axis([0 4000 -20 0])
title('最大熵法Order=20滤波后的信号功率谱');grid on
%%
%%功率谱密度
%PSD_WELCH方法
Nfft=512;n=0:sigLength-1;t=n/Fs; %数据长度、时间序列
window=hanning(256); %选用的窗口
noverlap=128; %分段序列重叠的采样点数（长度）
dflag='none'; %不做趋势处理
[Pxx,Pxxc,f]=psd(x,Nfft,Fs,window,noverlap,; %功率谱估计,并以的置信度给出置信区间，无返回值是绘制出置信区间
figure;subplot(211);
plot(f,10*log10(Pxx)); %绘制功率谱
xlabel('频率/Hz');ylabel('功率谱/dB');
grid on;title('PSD—Welch方法的原始信号功率谱')
subplot(212)
[Pxx,Pxxc,f]=psd(y0,Nfft,Fs,window,noverlap,; %功率谱估计,并以的置信度给出置信区间，无返回值是绘制出置信区间
plot(f,10*log10(Pxx)); %绘制功率谱
xlabel('频率/Hz');ylabel('功率谱/dB');axis([0 4000 -30 0])
grid on;title('PSD—Welch方法的滤波后的信号功率谱')
%%
%用多窗口法(MTM)
n=0:sigLength-1;t=n/Fs; %数据长度、分段数据长度，时间序列
[Pxx1,f]=pmtm(x,2,Nfft,Fs); %用多窗口法(NW=4)估计功率谱
figure;subplot(2,1,1),plot(f,10*log10(Pxx1)); %绘制功率谱
功率谱密度
振动台在使用中经常运用的公式 1、 求推力（F ）的公式 F=（m 0+m 1+m 2+ ……）A …………………………公式（1） 式中：F —推力（激振力）（N ） m 0—振动台运动部分有效质量（kg ） m 1—辅助台面质量（kg ） m 2—试件（包括夹具、安装螺钉）质量（kg ） A — 试验加速度（m/s 2） 2、 加速度（A ）、速度（V ）、位移（D ）三个振动参数的互换运算公式 2.1 A=ωv ……………………………………………………公式（2） 式中：A —试验加速度（m/s 2） V —试验速度（m/s ） ω=2πf （角速度） 其中f 为试验频率（Hz ） 2.2 V=ωD ×10-3 ………………………………………………公式（3） 式中：V 和ω与“2.1”中同义 D —位移（mm 0-p ）单峰值 2.3 A=ω2 D ×10-3 ………………………………………………公式（4） 式中：A 、D 和ω与“2.1”，“2.2”中同义 公式（4）亦可简化为： A= D f ?250 2 式中：A 和D 与“2.3”中同义，但A 的单位为g 1g=9.8m/s 2 所以： A ≈D f ?25 2 ,这时A 的单位为m/s 2 定振级扫频试验平滑交越点频率的计算公式 3.1 加速度与速度平滑交越点频率的计算公式 f A-V = V A 28.6 ………………………………………公式（5） 式中：f A-V —加速度与速度平滑交越点频率（Hz ）（A 和V 与前面同义）。
3.2 速度与位移平滑交越点频率的计算公式 D V f D V 28.6103?=- …………………………………公式（6） 式中：D V f -—加速度与速度平滑交越点频率（Hz ）（V 和D 与前面同义）。 3.3 加速度与位移平滑交越点频率的计算公式 f A-D =D A ??2 3 )2(10π ……………………………………公式（7） 式中：f A-D — 加速度与位移平滑交越点频率（Hz ），（A 和D 与前面同义）。 根据“3.3”，公式（7）亦可简化为： f A-D ≈5× D A A 的单位是m/s 2 4、 扫描时间和扫描速率的计算公式 4.1 线性扫描比较简单： S 1= 1 1 V f f H - ……………………………………公式（8） 式中： S1—扫描时间（s 或min ） f H -f L —扫描宽带，其中f H 为上限频率，f L 为下限频率（Hz ） V 1—扫描速率（Hz/min 或Hz/s ） 4.2 对数扫频： 4.2.1 倍频程的计算公式 n=2Lg f f Lg L H ……………………………………公式（9） 式中：n —倍频程（oct ） f H —上限频率（Hz ） f L —下限频率（Hz ） 4.2.2 扫描速率计算公式 R= T Lg f f Lg L H 2/ ……………………………公式（10） 式中：R —扫描速率（oct/min 或）
功率及功率谱计算
功率谱定义 从确定性信号功率计算开始 ()()221 11lim lim 222T T T T T P x t dt X d T T ωωπ∞--∞→∞→∞==?? ()()21lim 2T T S X T ωω→∞= S(w)为功率谱密度，简称功率谱 则 ()12P S d ωωπ+∞-∞= ? 随机信号的功率谱密度 （1）样本功率谱与功率谱密度 ()()21,lim ,2X T T S X T ωξωξ→∞= 针对一个具体的样本而言，其是一个确定性的信号 (2) 随机信号的平均功率及平均功率谱密度 ()X X P E P ξ=???? 需要对具体的样本取概率均值才能计算出功率 ()()()21,lim ,2X X T T S E S E X T ωωξωξ→∞??==?????? 故功率谱密度是对所有概率取期望的反应。 （3）自相关函数与功率谱密度 ()()R S τω? （4）信号的自相关函数计算 分为确定信号和随机信号 确定信号 02002*0 1()lim ()()T T x T R x t x t dt T ττ-→∞=-? 周期信号 0202*0 1()()()T T x R x t x t dt T ττ-=-? 随机信号 *()[()()]x R E x t x t ττ=- 2 功率计算 （1）根据定义来计算
（2）周期信号如何计算 0cos()A t ω的计算 200()()1()[]2 A A s d T πσωωπσωωωω+∞-∞-++==?不好算因此放弃，但是应该可以类推得出结论 (3)自相关函数计算 0cos()A t ω的计算 /2 200/2 /222000/2201()cos()cos(())cos()cos(2)1[]2 cos()2 T T T T r A t t d T A A t d T A τωωτωωτωωτωωτ+-+-=-+-==?? 所以其功率谱为 200()2 A πσωωσωω（-）+（+） 0j t Ae ω的计算 0000/2()2/2 /22/2 21()1T j t j t T T j T j r A e e dt T A e dt T A e ωωτωτωτ τ+---+-===?? 总结：因此周期函数，首先转换成傅里叶级数，然后再通过自相关函数的定义计算自相关函数，得到其功率谱密度。
功率谱密度
t=0:0.0001:0.1; %时间间隔为0.0001，说明采样频率为10000Hz x=square(2*pi*1000*t); %产生基频为1000Hz的方波信号 n=randn(size(t)); %白噪声 f=x+n; %在信号中加入白噪声 figure(1); subplot(2,1,1); plot(f); %画出原始信号的波形图 ylabel('幅值(V)'); xlabel('时间(s)'); title('原始信号'); y=fft(f,1000); %对原始信号进行离散傅里叶变换，参加DFT采样点的个数为1000 subplot(2,1,2); m=abs(y); f1=(0:length(y)/2-1)'*10000/length(y);%计算变换后不同点对应的幅值plot(f1,m(1:length(y)/2)); ylabel('幅值的模'); xlabel('时间(s)'); title('原始信号傅里叶变换'); %用周期图法估计功率谱密度 p=y.*conj(y)/1000; %计算功率谱密度 ff=10000*(0:499)/1000; %计算变换后不同点对应的频率值 figure(2); plot(ff,p(1:500)); ylabel('幅值'); xlabel('频率(Hz)'); title('功率谱密度（周期图法）'); 功率谱估计在现代信号处理中是一个很重要的课题，涉及的问题很多。在这里，结合matlab，我做一个粗略介绍。功率谱估计可以分为经典谱估计方法与现代谱估计方法。经典谱估计中最简单的就是周期图法，又分为直接法与间接法。直接法先取N点数据的傅里叶变换（即频谱），然后取频谱与其共轭的乘积，就得到功率谱的估计；间接法先计
功率谱估计真实验
功率谱估计仿真实验 选题条件：对于给定的一个信号()()()t t f t f t x ?ππ++=212sin 2)2sin(，其中1f =50Hz ， 2f =100Hz ，()t ?为白噪声，采样频率Fs 为1000Hz ，对其进行功率谱估 计。 仿真目标：采用多种方法对该指定信号进行功率谱估计，计算其功率谱密度，比较 各种估计方法的优劣。 设计思路：本仿真实验采用经典谱估计中的周期图法对给定信号进行谱估计。但是 由于其自身的缺陷，使得频率分辨率较低。为了不断满足需要，找到恰 当的估计法，实验使依次使用了周期图法的改进型方法如分段周期图法、 窗函数法以及修正的周期图法进行功率谱估计，对四种方法得出的谱估 计波形进行比较分析，得出估计效果最好的基于周期图法的谱估计方法。 仿真指标：频率分辨率、估计量的方差、频谱光滑度 平台说明：本实验采用MATLAB7.0仿真软件，基于WINDOWS-XP 系统。Matlab 是 一个集数值分析、矩阵运算、信号处理和图形显示于一体的工程分析处理软件。它提供的部分算法函数为功率谱估计提供了一条可行的方便途径，如PSD 和CSD 可以自动实现Welch 法估计，而不需要自己编程。但是较为有限，大部分需要自己编写相应的M 文件来实现。 实现方法： 一、周期图法 周期图法是直接将信号的采样数据()n x 进行傅立叶变换求功率谱密度估计。假设有限长随机信号序列()n x ，将它的功率谱按定义写出如下： ()()??? ?????+=∑-=-∞→2121lim N N n n j N j xx e n x N E e P ωω 如果忽略上式中求统计平均的运算，观测数据为：()n x 10-≤≤N n ，便得到了周期图法的定义： ()()2 10 ^ 1n j N n j xx e n x N e P ωω--=∑=， 式中的绝对值符号内的部分可以用FFT 计算，这样就可得到周期图法的计算框图如下所示： () ω j xx e ^ 图1 周期图法计算功率谱框图 采用周期图法时，可以分取不同的信号长度256、512和1024，分别进行功率谱
matlab实现功率谱密度分析psd
matlab实现功率谱密度分析psd及详细解说 功率谱密度幅值的具体含义？？ 求信号功率谱时候用下面的不同方法，功率谱密度的幅值大小相差很大！ 我的问题是，计算具体信号时，到底应该以什么准则决定该选用什么方法啊？ 功率谱密度的幅植的具体意义是什么？？下面是一些不同方法计算同一信号的matlab 程序！欢迎大家给点建议！ 直接法： 直接法又称周期图法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。 Matlab代码示例： 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)); 间接法： 间接法先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，便得到x(n)的功率谱估计。 Matlab代码示例： 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);
功率谱和功率谱密度的区别
谱让人联想到的Fourier变换，是一个时间平均（time average）概念，对能量就是能量谱，对功率就是功率谱。 功率谱的概念是针对功率有限信号的，所表现的是单位频带内信号功率随频率的变化情况。保留了频谱的幅度信息，但是丢掉了相位信息，所以频谱不同的信号其功率谱是可能相同的。 有两点需要注意： 1. 功率谱是随机过程的统计平均概念，平稳随机过程的功率谱是一个确定函数；而频谱是随机过程样本的Fourier变换，对于一个随机过程而言，频谱也是一个“随机过程”。（随机的频域序列） 2. 功率概念和幅度概念的差别。此外，只能对宽平稳的各态历经的二阶矩过程谈功率谱，其存在性取决于二阶矩是否存在并且二阶矩的Fourier变换收敛；而频谱的存在性仅仅取决于该随机过程的该样本的Fourier变换是否收敛。 频谱分析： 对动态信号在频率域内进行分析，分析的结果是以频率为坐标的各种物理量的谱线和曲线，可得到各种幅值以频率为变量的频谱函数F(ω)。频谱分析中可求得幅值谱、相位谱、功率谱和各种谱密度等等。频谱分析过程较为复杂，它是以傅里叶级数和傅里叶积分为基础的。 功率谱密度： 功率谱密度（PSD），它定义了信号或者时间序列的功率如何随频率分布。这里功率可能是实际物理上的功率，或者更经常便于表示抽象的信号被定义为信号数值的平方，也就是当信号的负载为1欧姆(ohm)时的实际功率。
由于平均值不为零的信号不是平方可积的，所以在这种情况下就没有傅里叶变换。维纳-辛钦定理（Wiener-Khinchin theorem）提供了一个简单的替换方法，如果信号可以看作是平稳随机过程，那么功率谱密度就是信号自相关函数的傅里叶变换。 信号的功率谱密度当且仅当信号是广义的平稳过程的时候才存在。如果信号不是平稳过程，那么自相关函数一定是两个变量的函数，这样就不存在功率谱密度，但是可以使用类似的技术估计时变谱密度。 随机信号是时域无限信号，不具备可积分条件，因此不能直接进行傅氏变换。一般用具有统计特性的功率谱来作为谱分析的依据。 功率谱与自相关函数是一个傅氏变换对。 功率谱具有单位频率的平均功率量纲。所以标准叫法是功率谱密度。从名字分解来看就是说，观察对象是功率，观察域是谱域。 通过功率谱密度函数，可以看出随机信号的能量随着频率的分布情况。像白噪声就是平行于一条直线。 一般我们讲的功率谱密度都是针对平稳随机过程的，由于平稳随机过程的样本函数一般不是绝对可积的，因此不能直接对它进行傅立叶分析。可以有三种办法来重新定义谱密度，来克服上述困难。 1. 用相关函数的傅立叶变换来定义谱密度； 2. 用随机过程的有限时间傅立叶变换来定义谱密度； 3. 用平稳随机过程的谱分解来定义谱密度。 三种定义方式对应于不同的用处，首先第一种方式前提是平稳随机过程不包含周
FFT在功率谱密度计算中的应用
F F T在功率谱密度计算 中的应用 集团标准化工作小组 #Q8QGGQT-GX8G08Q8-GNQGJ8-MHHGN#
FFT在功率谱密度计算中的应用 一、FFT算法理论依据和编程思想 FFT算法的基本思想： 考察DFT与IDFT的运算发现，利用以下两个特性可减少运算量： Ⅰ)系数是一个周期函数，它的周期性和对称性可利用来改进运算，提高计算效率。如: 因此 利用这些周期性和对称性，DFT运算中有些项可合并； Ⅱ）利用W N nk的周期性和对称性，把长度为N点的大点数的DFT运算分解为若干个小点数的DFT。因为DFT的计算量正比于N2，N小计算量也就小。 FFT算法正是基于这样的基本思想发展起来的。它有多种形式，下面是按时间抽取的FFT（N点DFT运算的分解） 先从一个特殊情况开始，假定N是2的整数次方，N=2M，M：正整数 1.将N点的DFT分解为两个N/2点的DFT： 首先将序列x（n）分解为两组，一组为偶数项，一组为奇数项 r=0,1，…，N/2-1 将DFT运算也相应分为两组： 其中X 1（k）和X 2 （k）分别是x 1 （r）和x 2 （r）的N/2点DFT。 可见，一个N点的DFT可以分解为两个N/2点的DFT，这两个N/2点的DFT再按照上面 （1）式合成为一个N点DFT，注意到，X 1（k），X 2 （k）有N/2个点，即k=0，1，…， N/2-1，由（1）式得到X（k）只有N/2点，而实际上X（k）有N个点，即k=0， 1，…，N-1，要用X 1（k），X 2 （k）表示全部X（K）值，还必须应用系数w的周期性和 对称性。 (k)的(N/2)～N-1点表示： 由X(k)= X 1（k）+w k N X 2 （k）, k=0,1,2,…,N/2-1
功率谱图应用
1.基本方法 周期图法是直接将信号的采样数据x(n)进行Fourier变换求取功率谱密度估计的方法。假定有限长随机信号序列为x(n)。它的Fourier变换和功率谱密度估计存在下面的关系： 式中，N为随机信号序列x(n)的长度。在离散的频率点f=kΔf,有： 其中，FFT[x(n)]为对序列x(n)的Fourier变换，由于FFT[x(n)]的周期为N，求得的功率谱估计以N为周期，因此这种方法称为周期图法。下面用例子说明如何采用这种方法进行功率谱 用有限长样本序列的Fourier变换来表示随机序列的功率谱，只是一种估计或近似，不可避免存在误差。为了减少误差，使功率谱估计更加平滑，可采用分段平均周期图法（Bartlett法）、加窗平均周期图法（Welch 法）等方法加以改进。 2. 分段平均周期图法（Bartlett法） 将信号序列x(n)，n=0,1,…,N-1，分成互不重叠的P个小段，每小段由m个采样值，则P*m=N。对每个小段信号序列进行功率谱估计，然后再取平均作为整个序列x(n)的功率谱估计。 平均周期图法还可以对信号x(n)进行重叠分段，如按2:1重叠分段，即前一段信号和后一段信号有一半是重叠的。对每一小段信号序列进行功率谱估计，然后再取平均值作为整个序列x(n)的功率谱估计。这两种方法都称为平均周期图法，一般后者比前者好。程序运行结果为图9-5，上图采用不重叠分段法的功率谱估计，下图为2:1重叠分段的功率谱估计，可见后者估计曲线较为平滑。与上例比较，平均周期图法功率谱估计具有明显效果（涨落曲线靠近0dB）。 3.加窗平均周期图法 加窗平均周期图法是对分段平均周期图法的改进。在信号序列x(n)分段后，用非矩形窗口对每一小段信号序列进行预处理，再采用前述分段平均周期图法进行整个信号序列x(n)的功率谱估计。由窗函数的基本知识（第7章）可知，采用合适的非矩形窗口对信号进行处理可减小“频谱泄露”，同时可增加频峰的宽度，从而提高频谱分辨率。 其中上图采用无重叠数据分段的加窗平均周期图法进行功率谱估计，而下图采用重叠数据分段的加窗平均周期图法进行功率谱估计，显然后者是更佳的，信号谱峰加宽，而噪声谱均在0dB附近，更为平坦（注意采用无重叠数据分段噪声的最大的下降分贝数大于5dB，而重叠数据分段周期图法噪声的最大下降分贝数小于5dB）。 4. Welch法估计及其MATLAB函数 Welch功率谱密度就是用改进的平均周期图法来求取随机信号的功率谱密度估计的。Welch 法采用信号重叠分段、加窗函数和FFT算法等计算一个信号序列的自功率谱估计(PSD如上例中的下半部分的求法)和两个信号序列的互功率谱估计（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为分段序列重叠的采样
功率谱密度 的估计
功率谱密度的估计 原始波=余弦波+白噪声 这个实验采用了两个输入，一个是白噪声，一个是有用信号和噪声信号作为输入时，他们的功率谱密度的仿真图像，并将他们进行对比。 平稳随机信号的功率谱密度（PSD ）是相关序列的离散傅里叶变换： ()()jw m XX x P w r m e ∞ --∞=∑ 采用间接法计算噪声信号的功率谱。 间接法，又称自相关法或者BT 法，在1985年由布莱克曼与图基首先开拓。间接法的理论基础是维纳-辛钦定理。他是由N 个观察值x(0),x(1),……,x(N-1),估计出自相关函数R （m ）,然后再求R （m ）的傅里叶变换作为功率谱密度的估计。 ()(),||1M jw jw m N m M S e R m e M N -=-=<=-∑ clear all; randn('state',0) NFFT=1024; %采样点数 Fs=1000; %取样频率（单位为Hz ） t=0:1/Fs:.2;
y1=cos(t*20*pi); %余弦序列 figure(1) plot(t,y1); ylabel('余弦序列'); grid on; %余弦序列的图像： %白噪声 m=(0:NFFT-1)/Fs; y=0.1*randn(size(m)); %产生高斯白噪声。 figure(2); plot(m,y); title('白噪声波形'); grid on;
%白噪声的自相关函数 [cory,lags]=xcorr(y,200,'unbiased'); %计算白噪声的自相关函数 figure(3) plot(lags,cory); %自相关函数(无偏差的)，其中，cory为要求的自相关函数，lag为自相关函数的长度。 title('白噪声相关函数'); grid on;
功率谱密度
功率谱密度 不同形式的数字基带信号具有不同的频谱结构，分析数字基带信号的频谱特性，以便合理地设计数字基带信号，使得消息代码变换为适合于给定信道传输特性的结构，是数字基带传输必须考虑的问题。 在通信中，除特殊情况（如测试信号）外，数字基带信号通常都是随机脉冲序列。因为，如果在数字通信系统中所传输的数字序列是确知的，则消息就不携带任何信息，通信也就失去了意义。故我们面临的是一个随机序列的谱分析问题。 考察一个二进制随机脉冲序列。设脉冲、分别表示二进制码“0”和“1”, 为 码元的间隔，在任一码元时间内，和出现的概率分别为p和1-p。 则随机脉冲序列x(t)可表示成: 其中 研究由上面二式所确定的随机脉冲序列的功率谱密度，要用到概率论与随机过程的有关知识。可以证明，随机脉冲序列x(t)的双边功率谱公式（1）： 其中、分别为、的傅氏变换，。 可以得出如下结论： （1）随机脉冲序列功率谱包括两部分：连续谱（第一项）和离散谱（第二项）。对于连续谱而言，由于代表数字信息的及不能完全相同，故，因此，连 续谱总是存在；而对于离散谱而言，则在一些情况下不存在，如及是双极性的脉冲，且出现概率相同时。 （2）当、、p及给定后，随机脉冲序列功率谱就确定了。 上式的结果是非常有意义的，它一方面能使我们了解随机脉冲序列频谱的特点，以及如何去具体地计算它的功率谱密度；另一方面根据它的离散谱是否存在这一特点，将使我们明确能否从脉冲序列中直接提取离散分量，以及采取怎样的方法可以从基带脉冲序列中获得所需的离散分量。这一点，在研究位同步、载波同步等问题时，将是十分重要的；再一方面，根据它的连续谱可以确定序列的带宽（通常以谱的第一个零点作为序列的带宽）。 下面，以矩形脉冲构成的基带信号为例，通过几个有代表性的特例对功率谱密度公式的应用及意义做进一步的说明，其结果对后续问题的研究具有实用意义。
滤波与功率谱估计
清华大学 《数字信号处理》期末作业 2013 年 1 月
第一题掌握去噪的方法 1.1 题目描述 MATLAB 中的数据文件noisdopp 含有噪声，该数据的抽样频率未知。调出该数据，用你学过的滤波方法和奇异值分解的方法对其去噪。要求：1．尽可能多地去除噪声，而又不损害原信号； 2．给出你去噪的原理与方法；给出说明去噪效果的方法或指标； 3．形成报告时应包含上述内容及必要的图形，并附上原程序。 1.2 信号特性分析 MATLAB所给noisdopp信号极其频域特征如图1.1、图1.2。 图1.1含有噪声的noisdopp信号
图1.2 noisdopp 信号频域特性 其中横坐标f 采用归一化频率，即未知抽样频率Fs 对应2（与滤波器设计时参数一致）。信号基本特性是一个幅值和频率逐渐增加的正弦信号叠加噪声，噪声为均匀的近似白噪声，没有周期等特点。 因为噪声信号能量在全频带均匀分布，滤波器截止频率过低则信号损失大，过高则噪声抑制小，认为频谱中含有毛刺较多的部分即为信噪比较小的部分，滤除这部分可以达到较好的滤波效果。 先给定去噪效果的评定指标。信号开始阶段频率较高（如图1.3，红圈为信号值），一周期内采样点4~5个，即信号归一化频率达到0.4~0.5（Fs=2），难以从频域将这部分信号同噪声分离，滤波后信号损失较大，故对前128点用信噪比考察其滤波效果，定义： 2 2 () 10lg (()())k k x k SNR y k x k =-∑∑ 其中，()x k 为原nosidopp 信号，()y k 为滤波后信号。SNR 越大表示滤除部分能力越小，可以反映滤波后信号对原信号的跟踪能力，对前128点主要考察SNR ，越大滤波器性能越好。
功率谱密度
功率谱密度谱是一种概率统计方法，是对随机变量均方值的量度。一般用于随机振动分析，连续瞬态响应只能通过概率分布函数进行描述，即出现某水平响应所对应的概率。 功率谱密度是结构在随机动态载荷激励下响应的统计结果，是一条功率谱密度值—频率值的关系曲线，其中功率谱密度可以是位移功率谱密度、速度功率谱密度、加速度功率谱密度、力功率谱密度等形式。数学上，功率谱密度值—频率值的关系曲线下的面积就是方差，即响应标准偏差的平方值。 谱是个很不严格的东西，常常指信号的Fourier变换，是一个时间平均（time average）概念功率谱的概念是针对功率有限信号的(能量有限信号可用能量谱分析)，所表现的是单位频带内信号功率随频率的变换情况。保留频谱的幅度信息，但是丢掉了相位信息，所以频谱不同的信号其功率谱是可能相同的。有两个重要区别：1。功率谱是随机过程的统计平均概念，平稳随机过程的功率谱是一个确定函数；而频谱是随机过程样本的Fourier变换，对于一个随机过程而言，频谱也是一个“随机过程”。（随机的频域序列）2。功率概念和幅度概念的差别。此外，只能对宽平稳的各态历经的二阶矩过程谈功率谱，其存在性取决于二阶局是否存在并且二阶矩的Fourier变换收敛；而频谱的存在性仅仅取决于该随机过程的该样本的Fourier变换是否收敛。热心网友回答提问者对于答案的评价：谢谢解答。 频谱分析（也称频率分析），是对动态信号在频率域内进行分析，分析的 结果是以频率为坐标的各种物理量的谱线和曲线，可得到各种幅值以频率为变 量的频谱函数F(ω)。频谱分析中可求得幅值谱、相位谱、功率谱和各种谱密 度等等。频谱分析过程较为复杂，它是以傅里叶级数和傅里叶积分为基础的。 功率谱是个什么概念？它有单位吗? 随机信号是时域无限信号，不具备可积分条件，因此不能直接进行傅氏变换。一般用具有统计特性的功率谱来作为谱分析的依据。功率谱与自相关函数是一个傅氏变换对。功率谱具有单位频率的平均功率量纲。所以标准叫法是功率谱密度。通过功率谱密度函数，可以看出随机信号的能量随着频率的分布情况。像白噪声就是平行于w轴，在w轴上方的一条直线。 功率谱密度，从名字分解来看就是说，观察对象是功率，观察域是谱域，通常指频域，密度，就是指观察对象在观察域上的分布情况。一般我们讲的功率谱密度都是针对平稳随机过程的，由于平稳随机过程的样本函数一般不是绝对可积的，因此不能直接对它进行傅立叶分析。可以有三种办法来重新定义谱密度，来克服上述困难。 一是用相关函数的傅立叶变换来定义谱密度；二是用随机过程的有限时间傅立叶变换来定义谱密度；三是用平稳随机过程的谱分解来定义谱密度。三种定义方式对应于不同的用处，首先第一种方式前提是平稳随机过程不包含周期分量并且均值为零，这样才能保证相关函数在时差趋向于无穷时衰减，所以lonelystar说的不全对，光靠相关函数解决不了许多问题，要求太严格了；对于第二种方式，虽然一个平稳随机过程在无限时间上不能进行傅立叶变换，但是对于有限区间，傅立叶变换总是存在的，可以先架构有限时间区间上的变换，在对时间区间取极限，这个定义方式就是当前快速傅立叶变换（FFT）估计谱密度的依据；第三种方式是根据维纳的广义谐和分析理论：Generalized harmonic analysis, Acta Math, 55(1930),117-258,利用傅立叶-斯蒂吉斯积分，对均方连续的零均值平稳随机过程进行重构，在依靠正交性来建立的。 另外，对于非平稳随机过程，也有三种谱密度建立方法，由于字数限制，功率谱密度的单位
FFT在功率谱密度计算中的应用
FFT在功率谱密度计算中的应用 一、FFT算法理论依据和编程思想 FFT算法的基本思想： 考察DFT与IDFT的运算发现，利用以下两个特性可减少运算量： Ⅰ)系数是一个周期函数，它的周期性和对称性可利用来改进运算，提高计算效率。如: 因此 利用这些周期性和对称性，DFT运算中有些项可合并； Ⅱ）利用W N nk的周期性和对称性，把长度为N点的大点数的DFT运算分解为若干个小点数的DFT。因为DFT的计算量正比于N2，N小计算量也就小。FFT算是基于这样的基本思想发展起来的。它有多种形式，下面是按时间抽取的FFT（N点DFT运算的分解）先从一个特殊情况开始，假定N是2的整数次方，N=2M，M：正整数 1.将N点的DFT分解为两个N/2点的DFT： 首先将序列x（n）分解为两组，一组为偶数项，一组为奇数项 r=0,1，…，N/2-1 将DFT运算也相应分为两组： 其中X1（k）和X2（k）分别是x1（r）和x2（r）的N/2点DFT。 可见，一个N点的DFT可以分解为两个N/2点的DFT，这两个N/2点的DFT再按照上面（1）式合成为一个N点DFT，注意到，X1（k），X2（k）有N/2个点，即k=0，1，…，N/2-1，由（1）式得到X（k）只有N/2点，而实际上X（k）有N个点，即k=0，1，…，N-1，要用X1（k），X2（k）表示全部X（K）值，还必须应用系数w的周期性和对称性。
2.X(k)的(N/2)～N-1点表示： 由X(k)= X1（k）+w k N X2（k）, k=0,1,2,…,N/2-1 （2a） 得： ， 因为 , 且 同样 。 考虑到W N k对称性：。 故（2b） （2a）式表示了X（k）前半部分k=0～N/2-1时的组成方式，（2b）式则表示了后半部分k=N/2～N-1时的组成方式。这两式所表示的运算过程可用一个称作蝶形的信号流图来表示。 3.蝶形信号流图：
噪声功率谱密度与方差之间的关系
关于matlab 中噪声功率谱密度与方差之间的关系的理解 1. 连续时间系统 高斯白噪声的定义为：如果一个噪声，它的幅度分布服从高斯分布，而它的功率谱密度又是均匀分布的，则称它为高斯白噪声。 故对于连续时间系统，理想的高斯白噪声的功率谱密度是一个常数，设为n0，而带宽是无限宽的，其功率为： 0*n ∞=∞ （1） 在n0不是为无穷小的情况下，理想的噪声功率Pn 是无限大的。 而实际当中，噪声带宽是有限宽的，只需要在我们所关心的频带范围内，噪声功率谱密度是个常数，则我们可认为其是高斯白噪声。设噪声单边功率谱密度为0n ，低通带宽为W ，则其噪声功率为： 0*2n n P W = （2） 如图1.1所示： o W -W 幅度 频率/HZ 0 2 n 图1.1 我们知道，高斯白噪声的分布为2 ~(0,)X N σ，则其功率为： 222()()()()n P E x D x E x D x σ==+== （3） 故对于低通系统有： 20/2 n W σ= （4） 而对于带通系统，如图1.2所示，有： 200*2*2n n P W n W σ=== （5）
W -W 幅度 频率/HZ 0 2 n 2. 离散时间系统 对于离散时间系统而言，带宽受到抽样速率fs 的限制。设WGN 一秒内抽取的一组数据样本为： 12[],,....fs x n x x x = 22([])0;([])([])E x n D x n E x n σ=== 2.1理论分析 由于时间为单个的离散点，故理想功率为0；但有下列定义：对于序列[]x n 的能量E 定义为序列各抽样值的平方和，则数据样本的能量为： 2221()*[()]*s f s s E x n f E x n f σ===∑ （6） 将功率定义为序列能量除以序列的时间，即 2*t s b E P f T σ==（单位：J/S ） （7） 式中，Tb 为序列时间，此处等于1S 。 如果功率单位采用W/symbol ，则有： 2/s t s P P f σ==（单位：J/symbol ） 2.2另一种理解 而实际当中，抽样点是一个时间段，认为1/s s T f =时间内的幅值就等于此抽样时刻的幅值，则单位抽样时间内的噪声能量为： 22***t s s s E E T f T σσ=== （6） 则噪声功率（单位：J/symbol ）为：
(完整word版)自己编写算法的功率谱密度的三种matlab实现方法
功率谱密度的三种matlab 实现方法 一：实验目的： (1)掌握三种算法的概念、应用及特点； (2)了解谱估计在信号分析中的作用； (3) 能够利用burg 法对信号作谱估计，对信号的特点加以分析。 二;实验内容： （1）简单说明三种方法的原理。 （2）用三种方法编写程序，在matlab 中实现。 （3）将计算结果表示成图形的形式，给出三种情况的功率谱图。 （4）比较三种方法的特性。 （5）写出自己的心得体会。 三：实验原理： 1.周期图法： 周期图法又称直接法。它是从随机信号x(n)中截取N 长的一段，把它视为能量有限x(n)真实功率谱)(jw x e S 的估计)(jw x e S 的抽样. 认为随机序列是广义平稳且各态遍历的，可以用其一个样本x(n)中的一段)(n x N 来估计该随机序列的功率谱。这当然必然带来误差。由于对)(n x N 采用DFT ，就默认)(n x N 在时域是周期的，以及)(k x N 在频域是周期的。这种方法把随机序列样本x(n)看成是截得一段)(n x N 的周期延拓，这也就是周期图法这个名字的来历。
2.相关法（间接法): 这种方法以相关函数为媒介来计算功率谱，所以又叫间接法。这种方法的具体步骤是： 第一步：从无限长随机序列x(n)中截取长度N 的有限长序列列 )(n x N 第二步：由N 长序列)(n x N 求（2M-1）点的自相关函数)(m R x ∧ 序列。 )()(1)(1 m n x n x N m R N n N N x += ∑-=∧ (2-1) 这里，m=-(M-1)…,-1,0,1…,M-1，M N ，)(m R x 是双边序列，但是由自相关函数的偶对称性式，只要求出m=0,。。。,M-1的傅里叶变换，另一半也就知道了。 第三步：由相关函数的傅式变换求功率谱。即 jwm M M m X jw x e m R e S ----=∧∧ ∑= )()(1) 1( 以上过程中经历了两次截断，一次是将x(n)截成N 长，称为加数据窗，一次是将x(n)截成（2M-1）长，称为加延迟窗。因此所得的功率谱仅是近似值，也叫谱估计，式中的)(jw x e S 代表估值。一般取M< 功率谱密度估计方法的MATLAB实现 在应用数学和物理学中，谱密度、功率谱密度和能量谱密度是一个用于信号的通用概念，它表示每赫兹的功率、每赫兹的能量这样的物理量纲。在物理学中，信号通常是波的形式，例如电磁波、随机振动或者声波。当波的频谱密度乘以一个适当的系数后将得到每单位频率波携带的功率，这被称为信号的功率谱密度（power spectral density, PSD）或者谱功率分布（spectral power distribution, SPD）。功率谱密度的单位通常用每赫兹的瓦特数（W/Hz）表示，或者使用波长而不是频率，即每纳米的瓦特数（W/nm）来表示。信号的功率谱密度当且仅当信号是广义的平稳过程的时候才存在。如果信号不是平稳过程，那么自相关函数一定是两个变量的函数，这样就不存在功率谱密度，但是可以使用类似的技术估计时变谱密度。信号功率谱的概念和应用是电子工程的基础，尤其是在电子通信系统中，例如无线电和微波通信、雷达以及相关系统。因此学习如何进行功率谱密度估计十分重要，借助于Matlab工具可以实现各种谱估计方法的模拟仿真并输出结果。下面对周期图法、修正周期图法、最大熵法、Levinson递推法和Burg法的功率谱密度估计方法进行程序设计及仿真并给出仿真结果。 以下程序运行平台：Matlab R2015a（8.5.0.197613） 一、周期图法谱估计程序 1、源程序 Fs=100000; %采样频率100kHz N=1024; %数据长度N=1024 n=0:N-1; t=n/Fs; xn=sin(2000*2*pi*t); %正弦波,f=2000Hz Y=awgn(xn,10); %加入信噪比为10db的高斯白噪声 subplot(2,1,1); plot(n,Y) title('信号') xlabel('时间');ylabel('幅度'); 功率谱密度幅值的具体含义？？ 求信号功率谱时候用下面的不同方法，功率谱密度的幅值大小相差很大！ 我的问题是，计算具体信号时，到底应该以什么准则决定该选用什么方法啊？ 功率谱密度的幅植的具体意义是什么？？下面是一些不同方法计算同一信号的matlab 程序！欢迎大家给点建议！ 一、直接法： 直接法又称周期图法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。 Matlab代码示例： 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)); 二、间接法： 间接法先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，便得到x(n)的功率谱估计。 Matlab代码示例： 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); 三、改进的直接法： 对于直接法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N太小，谱的分辨率又不好，因此需要改进。 FFT 在功率谱密度计算中的应用 一、FFT 算法理论依据和编程思想 FFT 算法的基本思想： 考察DFT 与IDFT 的运算发现，利用以下两个特性可减少运算量： Ⅰ)系数 是一个周期函数，它的周期性和对称性可利用来改进运算， 提高计算效率。 如: 因此 利用这些周期性和对称性，DFT 运算中有些项可合并； Ⅱ）利用W N nk 的周期性和对称性，把长度为N 点的大点数的DFT 运算分解为若干个小点数的DFT 。因为DFT 的计算量正比于N 2，N 小计算量也就小。 FFT 算法正是基于这样的基本思想发展起来的。它有多种形式，下面是按时间抽取的FFT （N 点DFT 运算的分解） 先从一个特殊情况开始，假定N 是2的整数次方，N=2M ，M ：正整数 1.将N 点的DFT 分解为两个N/2点的DFT ： 首先将序列x （n ）分解为两组，一组为偶数项，一组为奇数项 r=0,1，…，N/2-1 将DFT 运算也相应分为两组： 其中X 1（k ）和X 2（k ）分别是x 1（r ）和x 2（r ）的N/2点DFT 。 可见，一个N 点的DFT 可以分解为两个N/2点的DFT ，这两个N/2点的DFT 再按照上面（1）式合成为一个N 点DFT ，注意到，X 1（k ），X 2（k ）有N/2个点，即k=0，1，…， N/2-1，由（1）式得到X（k）只有N/2点，而实际上X（k）有N个点，即k=0，1，…， N-1，要用X 1（k），X 2 （k）表示全部X（K）值，还必须应用系数w的周期性和对称性。 2.X(k)的(N/2)～N-1点表示： 由X(k)= X 1（k）+w k N X 2 （k）, k=0,1,2,…,N/2-1 得： ， （2a） 因为 , 且 同样 。 考虑到W N k对称性：。 故（2b） （2a）式表示了X（k）前半部分k=0～N/2-1时的组成方式，（2b）式则表示了后半部分k=N/2～N-1时的组成方式。这两式所表示的运算过程可用一个称作蝶形的信号流图来表示。 2.2.3 功率谱密度 我们定义信号()t f 的能量（作用归一化处理）： 由电压()t f （或者电流()t f ）在Ω1电阻上消耗的能量： ?∞ ∞-=dt t f E )(2， （注释：22u R u i u E ==?=/） 积分值存在，信号的能量为有限值，称()t f 为能量信号。 对于能量无限大的信号（如周期性信号），我们考虑能量的时间平均值，这显然就是信号的平均功率。这种信号称作（平均）功率信号。 我们定义信号()t f 的平均功率，为电压()t f 在Ω1电阻上消耗的平均功率（简称功率）： ()?-∞→=22 21T T T dt t f T S lim 式中，T 是为求平均的时间区间。 为了更好地描述能量信号、功率信号，我们引入能量谱密度和功率谱密度概念。 能量谱密度、功率谱密度函数表示信号的能量、功率密度随频率变化的情况。 我们知道，非周期性信号的频谱宽度是无限的，然而，实际上信号的大部分功率是集中在某个有限的频谱宽度内。 通过研究功率谱密度，可以帮助了解信号的功率分布情况，确定信号的频带等。 对于能量信号()t f ，根据付里叶反变换有 ()()?∞+∞-ωωωπ =d e F t f t j 21 则信号的能量： ()()???∞∞-∞+∞-ω+∞∞-ωωπ ==dt d e F t f dt t f E t j ])[(21 2 ()()()()???∞+∞-∞+∞-∞+∞-ωωω-?ωπ =ω?ωπ=d F F d dt e t f F E t j *21 21 当()t f 为实信号时，)()(*ω=ωF F 。今后如无特别说明，都是指实信号， 电振动台在使用中经常运用的公式 1、 求推力（F ）的公式 F=（m 0+m 1+m 2+ ……）A …………………………公式（1） 式中：F —推力（激振力）（N ） m 0—振动台运动部分有效质量（kg ） m 1—辅助台面质量（kg ） m 2—试件（包括夹具、安装螺钉）质量（kg ） A — 试验加速度（m/s 2） 2、 加速度（A ）、速度（V ）、位移（D ）三个振动参数的互换运算公式 2.1 A=ωv ……………………………………………………公式（2） 式中：A —试验加速度（m/s 2） V —试验速度（m/s ） ω=2πf （角速度） 其中f 为试验频率（Hz ） 2.2 V=ωD ×10-3 ………………………………………………公式（3） 式中：V 和ω与“2.1”中同义 D —位移（mm 0-p ）单峰值 2.3 A=ω2 D ×10-3 ………………………………………………公式（4） 式中：A 、D 和ω与“2.1”，“2.2”中同义 公式（4）亦可简化为： A= D f ?250 2 式中：A 和D 与“2.3”中同义，但A 的单位为g 1g=9.8m/s 2 所以： A ≈D f ?25 2 ,这时A 的单位为m/s 2 定振级扫频试验平滑交越点频率的计算公式 3.1 加速度与速度平滑交越点频率的计算公式 f A-V = V A 28.6 ………………………………………公式（5） 式中：f A-V —加速度与速度平滑交越点频率（Hz ）（A 和V 与前面同义）。 3.2 速度与位移平滑交越点频率的计算公式 D V f D V 28.6103?=- …………………………………公式（6） 式中：D V f -—加速度与速度平滑交越点频率（Hz ）（V 和D 与前面同义）。 3.3 加速度与位移平滑交越点频率的计算公式 f A-D =D A ??2 3 )2(10π ……………………………………公式（7） 式中：f A-D — 加速度与位移平滑交越点频率（Hz ），（A 和D 与前面同义）。 根据“3.3”，公式（7）亦可简化为： f A-D ≈5× D A A 的单位是m/s 2 4、 扫描时间和扫描速率的计算公式 4.1 线性扫描比较简单： S 1= 1 1 V f f H - ……………………………………公式（8） 式中： S1—扫描时间（s 或min ） f H -f L —扫描宽带，其中f H 为上限频率，f L 为下限频率（Hz ） V 1—扫描速率（Hz/min 或Hz/s ） 4.2 对数扫频： 4.2.1 倍频程的计算公式 n=2Lg f f Lg L H ……………………………………公式（9） 式中：n —倍频程（oct ） f H —上限频率（Hz ） f L —下限频率（Hz ） 4.2.2 扫描速率计算公式 R= T Lg f f Lg L H 2/ ……………………………公式（10） 式中：R —扫描速率（oct/min 或）功率谱密度估计方法的MATLAB实现
计算功率谱密度
FFT在功率谱密度计算中的应用
通信技术概论信号的能量谱密度与功率谱密度
电振动台的振动功率谱密度计算