使用matlab计算时间序列的振幅谱、功率谱（PSD）以及幅度（dB）

基本区分

  • 振幅谱对应频率域的振幅 |F(f_1)|
  • 功率谱对应振幅的平方 psd_{f1} = [|F(f_1)|]^2
  • 幅度（dB）对应振幅平方的换算 := 10 lg(psd_{f1}) 【根据matlabdb2pow帮助文档】

使用FFT计算振幅谱和功率谱（Matlab）

fs = 100;

t = 0:1/fs:2-1/fs;
x = cos(2*pi*20*t) + randn(length(t));

N = length(x);
if(mod(N, 2) == 1)
N = N - 1;
end

fft_x = fft(x, N);
fft_x = fft_x(1:N/2+1);
ampsp_x = (1/N) * abs(fft_x);
psd_x = (1/N) * abs(fft_x).^2;

ampsp_x(2:end-1) = 2 * ampsp_x(2:end-1);
psd_x(2:end-1) = 2 * psd_x(2:end-1);

注意事项

  • FFT要求时间序列的长度为2的倍数
  • 由于时间序列f(x)是离散信号，对应的傅里叶变换F(f)为周期信号，虽然长度相同，但F(f)包含两个周期（负频率部分和正频率部分），因此，我们只需要正频率部分即可fft_x = fft_x(1:N/2+1);
  • 由于F(f)包含两个周期，而我们仅取了一个周期，故要恢复准确的振幅和功率谱，需要对振幅（功率谱）乘以2；但零频率和奈奎斯特频率（fs/2）不会出现两次，故这两个频率对应的值不必乘以2

matlab计算功率谱密度（PSD）

1、对于平稳信号，可以使用FFT计算PSD（代码如上）

2、对于非平稳信号，可以使用如短时傅里叶变换、Multitaper spectral estimate等工具计算【matlab中pmtm函数】

3、Multitaper spectral estimate对FFT的能量泄露进行了优化，能更准确稳定的求解功率谱（振幅谱直接开方即可amp=\sqrt{PSD}）

最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
平台声明：文章内容（如有图片或视频亦包括在内）由作者上传并发布，文章内容仅代表作者本人观点，简书系信息发布平台，仅提供信息存储服务。

推荐阅读更多精彩内容