浅谈快速傅里叶变化的Matlab与Python实现与区别

信号处理免不了要求频率、画频谱图，但Matlab的 fft() 函数与Python的 numpy.fft.fft() 与 scipy.fftpack.fft() 函数得到的是fft变化后的双边复数值，离画频谱图还有几句代码的距离。
基本原理不介绍了，下面直接懒人投喂，给出Matlab与Python的两个函数，直接调用即可画频谱图。

注：两种语言的fft算法是有区别的，最后细聊！

Matlab的fftlw函数

输入是信号序列、对应的时间序列、以及是否作图，输出可以得到单边归一化之后的频率与对应的振幅，通过输出可以直接画出常用的频谱图！

function [ F,M ] = fftlw( x,y,draw ) %FFTLW 快速傅里叶变化2021.10.26 %输入   x--时间 y--信号 draw--1为画频谱图，0为不画 %输出   F--频率 M--幅值   N=length(y);                       %采样点数 if(mod(N,2)>0)     N=N-1; end    Fs=(N-1)/(x(N)-x(1));              %采样频率 F=(N/2:N-1)*Fs/N-Fs/2 ;            %频率 y2=abs(fftshift(fft(y(1:N))));     %快速傅里叶变化 M=2*y2(N/2+1:N)/N;                 %归一化 M(1)=M(1)/2;                       %常量除以2  if draw==1                         %可视化     figure     plot(F,M)     xlabel('f/HZ')     ylabel('amplitude')     title('频谱图') end end

Python的fftlw函数

输入与matlab的略有点不同，分别是采样频率、信号序列、是否作图，输出与matlab的函数一致。

import numpy as np from scipy.fftpack import fft,ifft import matplotlib.pyplot as plt  def fftlw(Fs,y,draw):     '''     Parameters     ----------     Fs : 采样频率     y :  信号序列     draw :1为画频谱图，0为不画       Returns     -------     f : 频率     M : 幅值     '''     L=len(y)                              #采样点数     f = np.arange(int(L / 2)) * Fs / L    #频率     #M = np.abs(np.fft.fft(y))*2/L         #采用numpy.fft.fft()函数并归一化     M = np.abs((fft(y))) *2/L             #采用scipy.fftpack.fft()函数并归一化     M = M[0:int(L / 2)]                   #取单边谱     M[0]=M[0]/2                           #常量除以2          if draw==1:                           #可视化         plt.figure()         plt.rcParams['font.sans-serif']=['SimHei']         plt.rcParams['axes.unicode_minus'] = False         plt.plot(f,M)         plt.xlabel('f/HZ')         plt.ylabel('amplitude')         plt.title('频谱图')     return f,M

构造简单的信号对比两种语言fftlw效果

举个例子，构造如下信号验证所写函数的正确性：

⋅

sin

⁡

(

⋅

100

)

⋅

sin

⁡

(

⋅

200

)

y=3+t\cdot \sin (2\pi t\cdot 100)+3\cdot \sin (2\pi t\cdot 200)

$y = 3 + t \cdot sin (2 π t \cdot 100) + 3 \cdot sin (2 π t \cdot 200)$

其中，包括常数项，周期项和趋势项，理论上该信号的频率应该为0Hz、100Hz、200Hz（具体怎么算的补一补书知识）。在这里，我设置采样频率 fs=10000，产生10000个数据点，时域图如下：
浅谈快速傅里叶变化的Matlab与Python实现与区别 Matlab调用fftlw函数的主函数

fs=10000;                                     %采样频率 x=0:1/fs:(10000-1)/fs;                        %时间序列 y=sin(2*pi*x*100).*x+3*sin(2*pi*x*200)+3;     %信号序列 figure                                        %画时域图 plot(x,y) title('时域图') xlabel('t/s') ylabel('amplitude') [f,m]=fftlw(x,y,1);                           %快速傅里叶变化并画频谱图

得到的频谱图如下：发现0Hz、100Hz、200Hz处的幅值分别为3，0.5，3。0Hz与200Hz处的幅值完美对应时域中常数项与

(

⋅

200

)

sin (2\pi t\cdot 200)

$s i n (2 π t \cdot 200)$ 的系数；而100Hz项

(

⋅

200

)

sin (2\pi t\cdot 200)

$s i n (2 π t \cdot 200)$ 的系数不是常数而是

$t$ ，且时间是0-1s，该项傅里叶变化得到的是该段时间内的平均幅值，也就是0.5。
浅谈快速傅里叶变化的Matlab与Python实现与区别 Python调用fftlw函数的主函数
直接加在def fftlw()的后文调用他就行。

Fs=10000                #采用频率 L=10000                 #采样点数 t=np.arange(0,L/Fs,1/Fs)   #时间序列         y=np.sin(2*np.pi*t*100)*t+3*np.sin(2*np.pi*t*200)+3  #信号序列 f,M=fftlw(Fs,y,1)         #快速傅里叶变化并画频谱图