fft与psd的关系【傅里叶变化求功率谱】

2022-8-30到9-02学习笔记

  1. 1、查找python函数源码方法

函数名.file,返回函数在包中的位置,然后可以一级一级的追踪;
例子:

import matplotlib.pyplot as plt
print(plt.__file__)

返回:'C:\\Users\\Alien\\AppData\\Local\\Programs\\Python\\Python37\\lib\\site-packages\\matplotlib\\pyplot.py'

  1. 2、判断某些元素是否在某些元素中

借用numpy的cbook._check_in_list(),但是需要研究怎么用~

  • 补充_check_in_list的用法
# TODO 补充cbook._check_in_list怎么用
import matplotlib.cbook as cbook
  1. 3、判断是否为函数对象

callable()函数可以判断输入对象是否是函数
例子:

callable(123)
返回 False

def fun_help():
	return 1
callable(fun_help)
返回 True
  1. 4、plt.psd()的计算方法

(1)将数据按照NFFT的长度分块
(2)对每块数据进行滤波(窗函数默认hanning窗,窗系数大小为NFFT长度)
(3)每块数据*滤波器系数?????????
(3)对每块数据进行FFT处理
(4)用共轭的方法求功率
(5)功率/FS求出功率谱密度
(6)功率谱密度/窗函数功率
(7)用np.roll对结果进行翻转
(9)某一个频点的功率:用所有blk在该频点的平均值来替代

代码计算实例:

# 1、按照NFFT的大小把数据拆分成块,大小不够删掉尾部凑成整数倍
Nfft_len = 256  # psd的默认值是256
fs = 900
data_len = data_100.shape[0]
blk_num = data_len//Nfft_len
data_reshape = data_100[:blk_num*Nfft_len].reshape(blk_num,Nfft_len).T
# 2、对每块数据进行滤波(窗函数默认hanning窗,窗系数大小为NFFT长度)
window_coef = np.hanning(Nfft_len).reshape(-1,1)
# 3、每块数据*滤波器系数(因为滤波器长度和数据长度一样长,所以就直接相乘,相当于卷积)(滤波的目的是为了抗混叠,因为把数据直接拆成块,相当于加了矩形窗)
data_flt = data_reshape * window_coef
# data_flt = data_reshape * 1
# 4、对每块数据进行FFT处理
data_fft = np.fft.fft(data_flt,n=Nfft_len,axis=0)
# 5、用共轭的方法求功率
data_pow = abs(data_fft * np.conj(data_fft))   # 结果是功率,取abs是为了抛掉值为0的虚部
# 6、功率/FS求出功率谱密度
data_powfs = data_pow / fs
# 7、功率谱密度/窗函数功率
data_one_win = data_powfs / sum(abs(window_coef)**2)
# data_one_win = data_powfs / 1
# 8、用np.roll对结果进行翻转
data_fftshift = np.roll(data_one_win,-Nfft_len//2,axis=0)  # 垂直方向翻转:因为一个blk其实是一个频点周期
# 9、某一个频点的功率:用所有blk在该频点的平均值来替代
data_fft_psd = np.mean(data_fftshift,axis=1)
# 10、画图check
# fs_point = np.linspace(-fs/2,fs/2,Nfft_len)  # 错的
fs_point_temp = np.fft.fftfreq(Nfft_len, 1/fs)
fs_point = np.roll(fs_point_temp,-Nfft_len//2,axis=0)
plt.plot(fs_point,10*np.log10(data_fft_psd),'*-')  # 注意这里用的是real画图,是因为前面计算的结果是功率,所以这块虚部是0
plt.psd(data_100,Fs=fs)


# 总结:用代码简单归纳fft和psd的关系
data_fft0 = data_fft[:,0]
data_pow = abs(data_fft0)**2
data_pow_fs = data_pow/fs
data_pow_coef = data_pow_fs/sum(abs(window_coef)**2)
data_shift = np.roll(data_pow_coef,128)
plt.plot(fs_point, 10*np.log10(data_shift),'*-')

# 最后:psd的含义是功率谱,而fft的结果是每个频点幅值+相位,所以需要将fft的结果转换成功率然后/fs,
# 同时为了表示方便,将每个频点的功率转换成dB表示的形式

在这里插入图片描述

  1. 5、用noise产生信号有直流原因

根本原因是在生成noise信号时,用的是np.random.rand函数,它产生的数据是[0-1]均匀分布,所以数据的均值就是0.5,所以在0频附近会有直流,把生成的noise信号减去均值后,直流消失

# 生成底噪信号
noise_data = np.random.rand(8192)+np.random.rand(8192)*1j  # 生成的是[0-1]均匀分布,所以画出的频谱在0位置有尖
p1,f1 = plt.psd(noise_data,NFFT=1024,Fs=983.04)

在这里插入图片描述

data_mean = np.mean(noise_data)
data_diff_mean = noise_data - data_mean
plt.psd(data_diff_mean,NFFT=1024,Fs=983.04)

在这里插入图片描述

  1. 6、为什么原始带直流的noise信号用psd和scipy.singal.welch画出的结果不一样

根据原因就是两个函数对直流的默认处理不一样,体现在detrend的默认参数值,psd的默认是不处理,welch的默认是去直流

plt.subplots(3,1)
plt.subplot(311)
plt.psd(noise_data,NFFT=1024,Fs=983.04)
plt.subplot(312)
plt.psd(noise_data,NFFT=1024,Fs=983.04,detrend='constant')
plt.subplot(313)
f2,p2 = sig.welch(noise_data,983.40,window='hann',nperseg=1024,nfft=1024,return_onesided='False')  # 返回的是功率值
f2 = scipy.fft.fftshift(f2)
p2 = scipy.fft.fftshift(p2)
plt.plot(f2,10*np.log10(abs(p2)))  # 信号频响

在这里插入图片描述

  • 0
    点赞
  • 7
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
总结一下FFT和维纳辛钦定理PSD的问题-功率谱图.rar 早上在论坛上问了两个问题, 一个是关于FFT频谱时纵坐标的值问题 https://www.ilovematlab.cn/thread-27092-1-1.html 还有一个是用维纳辛钦定理PSD时出现的问题 https://www.ilovematlab.cn/thread-27133-1-1.html 经过达人们的指点,和自己的总结,获得一点心得,在这里与大家分享一下:) 1.FFT频谱 [CODE] Fs = 40; n = 0:1/Fs:159*1/Fs; x = sin sin; N = length; X = fftshift); Px1 = X.*conj/N; plot*Fs/N,Px1); grid on; axis title; 首先,fftshift的问题,以前上数字信号处理时,老师专门给提出了这个函数,但是我发现论坛里好多不太明白这个函数意义的,OO~,一般,fft得到的是频谱范围在【0-2*pi】范围内的频谱,以高频pi为中心,但是一般使用过程中,使用的频谱习惯以低频0为中心,fftshift的功能就是将频谱进行移位,使之在【-pi,pi】之间; 另外,纵坐标的问题,版主edifier2008提示说用/N的方法归一化,我试了一下,每次采样长度变大时,纵坐标的整体值都会变大,/N之后,值变为1之内了,但是并不是理论算法中得到的1. 图形如下: fft.jpg fft 2.维纳辛钦定理功率谱的问题 [CODE] Fs = 40; n = 0:1/Fs:159*1/Fs; x = sin sin; N = length; Rx = xcorr; Px2 = fftshift); plot*Fs/,abs); grid on; axis title; 图形如下: fftwei.jpg 程序中可以看出,也要使用fftshift对fft得到的频谱进行移位以得到以低频0为中心的频谱,另外,得到的功率谱纵轴值特别大,是不是也需要除以采样长度,我试了一下,仍然是很大,个人认为,在MATLAB中计算自相关函数以及计算FFT时,都没有对加和进行归一,将/N这一个系数可能都给省略掉了。 此外,我在很多教材里面看了不少里面的例题,都没有注意纵轴值的问题,我觉得在进行频谱分析,重点在于频率点,以及相近频率点的谱图是不是能够分辨出来,而对于各谱的大小,有个相对的比较即可。 不当之处,还望大家给与指正,:) :victory:
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值