用Matlab求离散信号的能量与功率怎么编程

关于谱能量，有这样一种解释，你可以试着去算一算

信号可以分成能量信号与功率信号，非周期能量信号具有能量谱密度，是傅立叶变换的平方，功率信号具有功率谱密度，其与自相关函数是一对傅立叶变换对，等于傅立叶变换的平方/区间长度。不能混淆。能量信号是没有功率谱的。

胡广书老师的书上找到这么一段话，“随机信号在时间上是无限的，在样本上也是无穷多，因此随机信号的能量是无限的，它应是功率信号。功率信号不满足付里叶变换的绝对可积的条件，因此其付里叶变换是不存在的。如确定性的正弦函数的付里叶变换是不存在，只有引入了冲激函数才求得其付里叶变换。因此，对随机信号的频谱分析，不再简单的是频谱，而是功率谱。”

对于确定性信号而言，里面存在能量信号，是没有功率谱密度的，也存在功率信号，是有功率谱密度的。所以信号的频谱与是否是确定性信号没有必然联系。

以下论点来源于研学论坛：

频谱是信号的傅立叶凳游变换。它描述了信号在各个频率上的分布大小。频谱的平方(当能量有限，平均功率为0时称为能量谱)描述了信号能量在各个频率上的分布大小。

计算过程中，都是通过样本数据的快速傅立叶变换来计算。但不同的是，信号的频谱是复数，包含幅频响应和相频响应，重复计算时的结果基本相同。 而随机信号的功率谱也可以对数据进行FFT，但必须计算模值衡粗嫌的平方，因为功率谱是实数。而且换一组样本后，计算的结果略有不同，因为随机信号的样本取值不同。要得到真实的功率谱必须进行多次平均，次数越多越好。

根据parseval定理，信号傅氏变换模平方被定义为能量谱，即单位频率范围内包含的信号能量。自然，能量跟功率有一个时间平均的关系，所以，能量谱密度在时间上平均就得到了功率谱。

matlab实现经典功率谱估计

fft做出来是频谱，psd做出来是功率谱；功率谱丢失了频谱的相位信息；频谱不同的信号其功率谱是可能相同的；功率谱是幅度取模后平方，结果是个实数

matlab中自功率谱密度直接用psd函数就可以求，按照matlab的说法，psd能实现Welch法估计，即相当于用改进的平均周期图法来求取咐手随机信号的功率谱密度估计。psd求出的结果应该更光滑吧。

1、直接法：

直接法又称周期图法，它是把随机序列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))

2、间接法：

间接法先由序列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)

3、改进的直接法：

对于直接法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N太小，谱的分辨率又不好，因此需要改进。

3.1、Bartlett法

Bartlett平均周期图的方法是将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

window=boxcar(length(n)) %矩形窗

noverlap=0 %数据无重叠

p=0.9 %置信概率

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

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

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)

这问题都懒得回答了，不同小波函数DWT后得到的小波系数个数会有差羡岩别，减半也是近似的，这是使用滤波器替代积分实现小波变换遇到的技术帆派搏细节问题，DB1是就是haar小波，其滤波器易于态祥设计，刚好可以完全精确减半，其它滴就难了，参见

---