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

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（具體怎么算的補一補書知識）。在這里，我設(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)#快速傅里葉變化并畫頻譜圖

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

采用實際的振動信號對比兩種語言fftlw效果

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

Matlab畫頻譜圖

Python畫頻譜圖

python的頻譜圖的幅值與原始數(shù)據(jù)量級差別較大，與matlab的頻譜圖也毫不相關(guān)，可能是底層傅里葉變換的算法不同所致，具體哪個正確還帶進一步考證?。?！

到此這篇關(guān)于Python與Matlab實現(xiàn)快速傅里葉變化的區(qū)別的文章就介紹到這了,更多相關(guān)Python 傅里葉變化內(nèi)容請搜索本站以前的文章或繼續(xù)瀏覽下面的相關(guān)文章希望大家以后多多支持本站！

版權(quán)聲明：本站文章來源標(biāo)注為YINGSOO的內(nèi)容版權(quán)均為本站所有，歡迎引用、轉(zhuǎn)載，請保持原文完整并注明來源及原文鏈接。禁止復(fù)制或仿造本網(wǎng)站，禁止在非maisonbaluchon.cn所屬的服務(wù)器上建立鏡像，否則將依法追究法律責(zé)任。本站部分內(nèi)容來源于網(wǎng)友推薦、互聯(lián)網(wǎng)收集整理而來，僅供學(xué)習(xí)參考，不代表本站立場，如有內(nèi)容涉嫌侵權(quán)，請聯(lián)系alex-e#qq.com處理。