快捷導(dǎo)航

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

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

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

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

Matlab的fftlw函數(shù)

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

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


N=length(y);                       %采樣點(diǎn)數(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的略有點(diǎn)不同，分別是采樣頻率、信號(hào)序列、是否作圖，輸出與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 :  信號(hào)序列
    draw :1為畫(huà)頻譜圖，0為不畫(huà) 

    Returns
    -------
    f : 頻率
    M : 幅值
    '''
    L=len(y)                              #采樣點(diǎn)數(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)造簡(jiǎn)單的信號(hào)對(duì)比兩種語(yǔ)言fftlw效果

舉個(gè)例子，構(gòu)造如下信號(hào)驗(yàn)證所寫(xiě)函數(shù)的正確性：

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

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

在這里插入圖片描述

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

fs=10000;                                     %采樣頻率
x=0:1/fs:(10000-1)/fs;                        %時(shí)間序列
y=sin(2*pi*x*100).*x+3*sin(2*pi*x*200)+3;     %信號(hào)序列
figure                                        %畫(huà)時(shí)域圖
plot(x,y)
title('時(shí)域圖')
xlabel('t/s')
ylabel('amplitude')
[f,m]=fftlw(x,y,1);                           %快速傅里葉變化并畫(huà)頻譜圖

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

在這里插入圖片描述

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

Fs=10000                #采用頻率
L=10000                 #采樣點(diǎn)數(shù)
t=np.arange(0,L/Fs,1/Fs)   #時(shí)間序列        
y=np.sin(2*np.pi*t*100)*t+3*np.sin(2*np.pi*t*200)+3  #信號(hào)序列
f,M=fftlw(Fs,y,1)         #快速傅里葉變化并畫(huà)頻譜圖

在這里插入圖片描述

圖和matlab的一模一樣！但是！

采用實(shí)際的振動(dòng)信號(hào)對(duì)比兩種語(yǔ)言fftlw效果

數(shù)據(jù)來(lái)源于西儲(chǔ)大學(xué)轉(zhuǎn)子實(shí)驗(yàn)臺(tái)振動(dòng)信號(hào)，采樣頻率為12000Hz，現(xiàn)取正常狀態(tài)下、轉(zhuǎn)速1796 rpm軸承振動(dòng)信號(hào)1000個(gè)點(diǎn)如下。粗略的觀察，有一個(gè)低頻信號(hào)大概周期為0.035 s，頻率大概 29Hz。

在這里插入圖片描述