Python下FFT分析图文频率分辨率与什么有关

1. FFT相关理论

1.1 离散傅里叶变换(DFT)

离散傅里叶变换(discrete Fourier transform) 傅里叶分析方法是信号分析的最基本方法,傅里叶变换是傅里叶分析的核心,通过它把信号从时间域变换到频率域,进而研究信号的频谱结构和变化规律。但是它的致命缺点是:计算量太大,时间复杂度太高,当采样点数太高的时候,计算缓慢,由此出现了DFT的快速实现,即下面的快速傅里叶变换FFT。

1.2 快速傅里叶变换(FFT)

计算量更小的离散傅里叶的一种实现方法。

1.3 归一化和取一半处理

比如有一个信号如下:
Y=A1+A2cos(2πω2+φ2)+A3cos(2πω3+φ3)+A4*cos(2πω4+φ4)
经过FFT之后,得到的“振幅图”中,
第一个峰值(频率位置)的模是A1的N倍,N为采样点,本例中为N=1400,此例中没有,因为信号没有常数项A1
第二个峰值(频率位置)的模是A2的N/2倍,N为采样点,
第三个峰值(频率位置)的模是A3的N/2倍,N为采样点,
第四个峰值(频率位置)的模是A4的N/2倍,N为采样点,
依次下去…
考虑到数量级较大,一般进行归一化处理,既然第一个峰值是A1的N倍,那么将每一个振幅值都除以N即可
FFT具有对称性,一般只需要用N的一半,前半部分即可。

2. FFT之后每个点的频率分辨率与采样频率(Fs)和采样点数(N)之间的关系

FFT频率分辨率=1/采样总时间=1/(N/Fs) = Fs/N(Hz)。
例如采样频率10k,采样点数1024,则FFT之后,频率分辨率=10k/1024=9.765625Hz。

3. 代码

import numpy as np
from scipy.fftpack import fft,ifft
import matplotlib.pyplot as plt
from matplotlib.pylab import mpl
 
mpl.rcParams['font.sans-serif'] = ['SimHei']   #显示中文
mpl.rcParams['axes.unicode_minus']=False       #显示负号

Fs = 1000                       #采样频率
N = 4096                        #采样的总点数
Time = N/Fs;                       #采样时间N/Fs


x = np.arange(0,Time,1.0/Fs)        #采样的x等差数列,0,Ts,2Ts,...,Time-Ts
half_x = x[range(int(N/2))]     #取一半区间

k = np.arange(N)
fx = k / Time;                  #x轴归一化到采样频率。FFT频率分辨率=1/采样总时间(Hz),即采样时间为0.2s,则频率分辨率5Hz。即采样时间为5s,则频率分辨率0.2Hz。
half_fx = fx[range(int(N/2))]
   
 
#设置需要采样的信号,频率分量有200,400和600
y=1+2*np.sin(2*np.pi*200*x) + 3*np.sin(2*np.pi*400*x)


 
fft_y=fft(y)                          #快速傅里叶变换
 
 
abs_y=np.abs(fft_y)                # 取复数的绝对值,即复数的模(双边频谱)
angle_y=np.angle(fft_y)            #取复数的角度
normalization_y=abs_y/N            #归一化处理(双边频谱)                              
normalization_half_y = normalization_y[range(int(N/2))]      #由于对称性,只取一半区间(单边频谱)
 
plt.subplot(231)
plt.plot(x,y)
plt.xlabel('t/s')
plt.ylabel('Amplitude')
plt.title('原始波形')
plt.tight_layout() #设置默认的间距
 
plt.subplot(232)
plt.plot(fx,fft_y,'black')
plt.xlabel('Freq/Hz')
plt.ylabel('Amplitude振幅')
plt.title('双边振幅谱(未求振幅绝对值)',fontsize=9,color='black') 
 
plt.subplot(233)
plt.plot(fx,abs_y,'r')
plt.xlabel('Freq/Hz')
plt.ylabel('Amplitude振幅')
plt.title('双边振幅谱(未归一化)',fontsize=9,color='red') 
 
plt.subplot(234)
plt.plot(fx,angle_y,'violet')
plt.xlabel('Freq/Hz')
plt.ylabel('phase相位')
plt.title('双边相位谱(未归一化)',fontsize=9,color='violet')

##plt.subplot(234)
##plt.plot(x[0:50],y[0:50],'violet')
##plt.xlabel('t/s')
##plt.ylabel('Amplitude')
##plt.title('原始波形(部分))',fontsize=9,color='violet')

plt.subplot(235)
plt.plot(fx,normalization_y,'g')
plt.xlabel('Freq/Hz')
plt.ylabel('Amplitude振幅')
plt.title('双边振幅谱(归一化)',fontsize=9,color='green')
 
plt.subplot(236)
plt.plot(half_fx,normalization_half_y,'blue')
plt.xlabel('Freq/Hz')
plt.ylabel('Amplitude振幅')
plt.title('单边振幅谱(归一化)',fontsize=9,color='blue')
 
plt.show()

4. 效果

在这里插入图片描述

MATLAB中使用FFT做频谱分析频率分辨率问题-频率分辨率.rar MATLAB中使用FFT做频谱分析频率分辨率问题 最近做FFT时,使用的采样频率和信号长度的取舍一直没有搞清楚,后来在论坛上发了一个贴子《总结一下使用FFT和维纳-辛钦定理求解PSD问题》(讨论见https://www.ilovematlab.cn/thread-27150-1-1.html,特别感谢会员songzy41,他的问题给了我很大启示),跟帖中给了我不少启示,并且让我对“频率分辨率”这个概念有了更深入的理解。再次一并感谢论坛的高手们。 频率分辨率,顾名思义,就是将信号中两个靠的很近的频谱分开的能力。 信号x长度为Ts,通过傅氏变换后得到X,其频率分辨率为Δf=1/T(Hz),若经过采样后,假设采样频率为fs=1/Ts,而进行频谱分析时要将这个无穷长的序列使用窗函数截断处理,假设使用矩形窗,我们知道,矩形窗的频谱为sinc函数,主瓣宽度可以定义为2*pi/M,M为窗宽,那么,时域相乘相当于频域卷积,频域内,这一窗函数能够分辨出的最近频率肯定不可能小于2*pi/M了,也就是如果数据长度不能满足2*pi/M<|w2-w1|(w2,w1为两个靠的很近的频率),那么在频谱分析时,频谱上将不能分辨出这两个谱,由于w2-w1=2*pi/fs=2*pi*Δf/fs也就是2*pi/M<2*piΔf/fs,得到Δf的限制为fs/M,这就是窗函数宽度的最小选择,就是说,根据Shannon采样定理确定了采样频率后,要根据靠的最近的谱峰来确定最小的采样长度,这样,所作出来的频谱才能分辨出那两个谱峰,也就是拥有了相应的频率分辨率。 几个例子: 考虑双正弦信号:x = sin sin;根据Shannon采样定理,采样频率要大于截止频率的两倍,这里选采样频率为80,那么,我们可以看到,Δf为0.2Hz,那么,最小的数据长度为0.2/80=400,但是对正弦信号的频谱分析经验告诉我们,在截断时截断时的数据要包含整周期,并且后面不宜补零以避免频谱泄露(这一点见胡广书《数字信号处理导论》,清华大学出版社),那么,我们要选择至少980个点,才能保含到一个整周期,另外,FFT的经验告诉我们作分析时最好选择2的整数次幂,我们选择靠的最近的1024点。分析结束。 [CODE] Fs = 80; n = 0:1/Fs:1023*1/Fs; x = sin sin; N = length; figure; X = fftshift); plot*Fs/N,abs*2/N); grid on; axis; 这是按照我们的分析进行的编程和图形 zheng.jpg 可以看出这两个谱峰很好的被分辨开来,9.8Hz不在谱线上,所以幅值不为1,以下是一些对比: [CODE] Fs = 80; n = 0:1/Fs:1023*1/Fs; x = sin sin; N = length; X = fftshift); figure; subplot plot*Fs/N,abs*2/N); grid on; axis; title; n = 0:1/Fs:979*1/Fs; x = sin sin; N = length; X = fftshift); subplot plot*Fs/N,abs*2/N); grid on; axis; title; n = 0:1/Fs:399*1/Fs; x = sin sin; N = length; X = fftshift); figure; subplot plot*Fs/N,abs*2/N); grid on; axis; title; Fs = 20; n = 0:1/Fs:1024*1/Fs; x = sin sin; N = length; X = fftshift); subplot plot*Fs/N,abs*2/N); grid on; axis; title; 结果如下: 1024.jpg 400.jpg 这是我在做FFT以及论坛中的问题时所得到的一点启发,不当之处还请大家指正。OO~ 频率分辨率.rar 为了方便大家,我将doc版报告和m文件一起上传,和帖子内容一样。OO~
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值