Python與Matlab實現(xiàn)快速傅里葉變化的區(qū)別

更新時間：2021年10月28日 15:22:49 作者：易~lw

信號處理免不了要求頻率、畫頻譜圖，但Matlab的fft()函數(shù)與Python的numpy.fft.fft()與scipy.fftpack.fft()函數(shù)得到的是fft變化后的雙邊復(fù)數(shù)值，離畫頻譜圖還有幾句代碼的距離?；驹聿唤榻B了，下面直接懶人投喂，給出Matlab與Python的兩個函數(shù)，直接調(diào)用即可畫頻譜圖

注：兩種語言的fft算法是有區(qū)別的，最后細(xì)聊！

Matlab的fftlw函數(shù)

輸入是信號序列、對應(yīng)的時間序列、以及是否作圖，輸出可以得到單邊歸一化之后的頻率與對應(yīng)的振幅，通過輸出可以直接畫出常用的頻譜圖！

function [ F,M ] = fftlw( x,y,draw )
%FFTLW 快速傅里葉變化2021.10.26
%輸入   x--時間 y--信號 draw--1為畫頻譜圖，0為不畫
%輸出   F--頻率 M--幅值


N=length(y);                       %采樣點數(shù)
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函數(shù)

輸入與matlab的略有點不同，分別是采樣頻率、信號序列、是否作圖，輸出與matlab的函數(shù)一致。

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)                              #采樣點數(shù)
    f = np.arange(int(L / 2)) * Fs / L    #頻率
    #M = np.abs(np.fft.fft(y))*2/L         #采用numpy.fft.fft()函數(shù)并歸一化
    M = np.abs((fft(y))) *2/L             #采用scipy.fftpack.fft()函數(shù)并歸一化
    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

構(gòu)造簡單的信號對比兩種語言fftlw效果

舉個例子，構(gòu)造如下信號驗證所寫函數(shù)的正確性：

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

其中，包括常數(shù)項，周期項和趨勢項，理論上該信號的頻率應(yīng)該為0Hz、100Hz、200Hz（具體怎么算的補(bǔ)一補(bǔ)書知識）。在這里，我設(shè)置采樣頻率 fs=10000，產(chǎn)生10000個數(shù)據(jù)點，時域圖如下：

在這里插入圖片描述

Matlab調(diào)用fftlw函數(shù)的主函數(shù)

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);                           %快速傅里葉變化并畫頻譜圖

得到的頻譜圖如下：發(fā)現(xiàn)0Hz、100Hz、200Hz處的幅值分別為3，0.5，3。0Hz與200Hz處的幅值完美對應(yīng)時域中常數(shù)項與s i n ( 2 π t ⋅ 200 ) 的系數(shù)；而100Hz項sin(2πt⋅200)的系數(shù)不是常數(shù)而是 t，且時間是0-1s，該項傅里葉變化得到的是該段時間內(nèi)的平均幅值，也就是0.5。

在這里插入圖片描述

Python調(diào)用fftlw函數(shù)的主函數(shù)
直接加在def fftlw()的后文調(diào)用他就行。

Fs=10000                #采用頻率
L=10000                 #采樣點數(shù)
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)         #快速傅里葉變化并畫頻譜圖

在這里插入圖片描述