matlab實現離散傅立葉變換及低通濾波
阿新 • • 發佈:2019-01-12
如圖感測器無濾波狀態下FZ資料為下列
匯入matlab使用工具箱分析圖如下:
將資料匯入matlab程式碼
clear;clc;close all load('data_nofliter') Fs=100; % 採集頻率 T=1/Fs; % 採集時間間隔 %訊號長度最好為偶數 N=length(data_z); % 採集訊號的長度 y = data_z; %將資料放入y向量 t=(0:1:N-1)*T; % 定義整個採集時間點 t=t'; % 轉置成列向量 % 繪製時域訊號 figure plot(t,y) xlabel('時間') ylabel('訊號值') title('時域訊號') % fft變換 Y=fft(y); % Y為fft變換結果,複數向量 Y_half=Y(1:N/2+1); % 只看變換結果的一半即可 A=abs(Y_half); % 複數的幅值(模) f=(0:1:N/2)*Fs/N; % 生成頻率範圍 f=f'; % 轉置成列向量 % 幅值修正 A_adj=zeros(N/2+1,1); A_adj(1)=A(1)/N; % 頻率為0的位置 A_adj(end)=A(end)/N; % 頻率為Fs/2的位置 A_adj(2:end-1)=2*A(2:end-1)/N; % 繪製頻率幅值圖 figure subplot(3,2,1) plot(f,A_adj) xlabel('頻率 (Hz)') ylabel('幅值 (修正後)') title('FFT變換幅值圖') grid on % 繪製頻譜相點陣圖 subplot(3,2,2) phase_angle=angle(Y_half); % angle函式的返回結果為弧度 phase_angle=rad2deg(phase_angle); plot(f,phase_angle) xlabel('頻率 (Hz)') ylabel('相位角 (degree)') title('FFT變換相點陣圖') grid on %%計算理想低通濾波器的頻率響應 fc= 15; H = zeros(length(f),1); H(f<fc)=1; HA = abs(H); Hangle = angle(H); Hangle = rad2deg(Hangle); %將弧度轉換為角度 subplot(3,2,3) plot(f,HA,'m') xlabel('頻率(Hz)') ylabel('H的幅值') grid on subplot(3,2,4) plot(f,Hangle,'m') xlabel('頻率(Hz)') ylabel('H的相位角') grid on %%計算輸出訊號的傅立葉變換 output_half = Y_half.*H; output_halfA=abs(output_half ); % 幅值修正 output_halfA_adj=zeros(N/2+1,1); output_halfA_adj(1)=output_halfA(1)/N; % 頻率為0的位置 output_halfA_adj(end)=output_halfA(end)/N; % 頻率為Fs/2的位置 output_halfA_adj(2:end-1)=2*output_halfA(2:end-1)/N; subplot(3,2,5) plot(f,output_halfA_adj,'r') xlabel('頻率 (Hz)') ylabel('Y的幅值') grid on output_halfangle=angle(output_half); output_halfangle=rad2deg(output_halfangle); % 將弧度轉換為角度,方便觀察 subplot(3,2,6) plot(f,output_halfangle,'r') xlabel('頻率 (Hz)') ylabel('Y的相位角 (degree)') grid on %% 輸出訊號的半譜補全成全譜 second_half=output_half(2:end-1); second_half=conj(second_half); second_half=flip(second_half); output=[output_half;second_half]; %% 對output做傅立葉反變換,得到時域輸出訊號y(t) OUTPUT=ifft(output); figure plot(t,OUTPUT,'r') xlabel('時間 /s') ylabel('訊號值') title('時域輸出訊號 y(t)') grid on
執行得到下圖
低通濾波的頻率為15HZ,資料抖動範圍有所縮小。