Python与Matlab实现快速傅里叶变化的区别(频谱频率)

注：两种语言的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效果

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

y=3+t⋅sin(2πt⋅100)+3⋅sin(2πt⋅200)

其中，包括常数项，周期项和趋势项，理论上该信号的频率应该为0Hz、100Hz、200Hz（具体怎么算的补一补书知识）。在这里，我设置采样频率 fs=10000，产生10000个数据点，时域图如下：

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处的幅值完美对应时域中常数项与s i n ( 2 π t ⋅ 200 ) 的系数；而100Hz项sin(2πt⋅200)的系数不是常数而是t，且时间是0-1s，该项傅里叶变化得到的是该段时间内的平均幅值，也就是0.5。

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)#快速傅里叶变化并画频谱图