DFT和FFT功率（波数）谱分析

2024-06-26 07:26| 来源: 网络整理| 查看: 265

参考帖子 https://www.mathworks.com/matlabcentral/answers/24965-fft2-function-in-matlab 一维傅里叶变换 Y(omega)=【f(t)exp(-iomega*t)*dt】【】表示负无穷到正无穷积分对时间变换得到的是频率谱，omega=2 * pi * f，为圆频率，f为频率，单位s-1，即赫兹，单位时间内的cycle数对空间变换得到的是波数谱，omega=2 * pi * k，k为波数，单位为m-1，单位空间内的cycle数编程时积分改为累加，即离散傅里叶变换（DFT）

matlab里有fft函数，采用快速傅里叶变换（FFT）方法，大大节省计算量 FFT原理为将2的n次方长度（缺少补零）的序列分为奇数列和偶数列，通过它们之间的关系，将序列不断一分为二，最后得到最短的序列再用DFT方法求解。两者得到的结果是一致的。

matlab 还有pwelch函数用于计算功率谱，其中采用加窗、分段（各段之间有重叠）平均、fft方法求解功率谱，相比直接采用fft方法，增加了自由度、可靠性、并减少了能量泄露和旁瓣效应（离散采样导致的，峰值功率减弱，并在两侧出现小的波动峰值），并且提供了置信度区间计算。

下面的程序分别采用DFT、FFT和pwelch方法求解同一时间序列功率谱，并进行了比较

%% DFT clear;clc;close all; Fs=1/86400; % sample frequency: 1 day n=0:1023; % number of points t=n/Fs; % time per1=4; % periodicity 1: 4 days per2=8; data=cos(2*pi/(86400*per1)*t)+2*cos(2*pi/(86400*per2)*t)+randn(size(t)); %% *dt?? figure plot(t/86400,data) xlabel('time(days)'); ylabel('amplitude'); title('data'); Nt=length(n); % number of points dt=1/Fs; % time increment f_N=Fs/2; % Nyquist frequency df=Fs/Nt; % Frequency increment freq=-f_N:df:f_N-df; % Frequency range (centered) data_frequency = zeros(size(data)); % Compute 1D Discrete Fourier Transform for j1 = 1 : Nt for j2 = 1 : Nt data_frequency(j1) = data_frequency(j1) + ... data(j2)*exp(-1i*(2*pi)*(freq(j1)*t(j2))); end end %此处 j1为对应omega的下标，j2为对应t的下标 % here we calculate centered frency range, convert it to one-sided frequency % range plot_pdata=10*log10(abs(data_frequency(Nt/2+1:end)).^2/Nt); plot_freq=freq(Nt/2+1:end)*86400; figure plot(plot_freq,plot_pdata); xlabel('frequency(cpd)');ylabel('Spectral /dB'); title('DFT'); %% FFT pdata=fft(data); % FFT calculate two-sided frency range, convert it to one-sided frequency % range figure plot_pdata=10*log10(abs(pdata(1:Nt/2).^2/Nt)); plot_freq=freq(1:Nt/2)*86400; plot(plot_freq,plot_pdata(1:Nt/2)); xlabel('freq');ylabel('Spectral /dB'); title('FFT'); %% function pwelch % window=boxcar(Nt); %矩形窗 % window=hamming(Nt); %海明窗 % window=blackman(Nt); %blackman窗 noverlap=[]; %数据重叠 window=200; range='onesided'; %频率间隔为[0 Fs/2]，只计算一半的频率 [pdata,fout]=pwelch(data,window,noverlap,Nt,Fs,range); plot_pdata=10*log10(pdata/Nt); figure X=fout*86400; plot(X,plot_pdata); xlabel('freq');ylabel('Spectral /dB'); title('pwelch');

二维傅里叶变换

%二维数据，列代表空间，行代表时间，二维傅里叶变换后得到频率-波数关系即频散关系。 %此处用随机数据分析，fftshift后，以0作为频率中心 %给出了DFT和fft2函数之间的比较，在精度允许的范围内，两者结果一致 Nx = 64; % Number of samples collected along first dimension Nt = 32; % Number of samples collected along second dimension dx = 1; % Distance increment (i.e., Spacing between each column) dt = .1; % Time increment (i.e., Spacing between each row) x = 0 : dx : (Nx-1)*dx; % distance t = 0 : dt : (Nt-1)*dt; % time data_spacetime = randn(Nt,Nx); % random 2D matrix Nyq_k = 1/(2*dx); % Nyquist of data in first dimension Nyq_f = 1/(2*dt); % Nyquist of data in second dimension dk = 1/(Nx*dx); % Wavenumber increment df = 1/(Nt*dt); % Frequency increment k = -Nyq_k : dk : Nyq_k-dk; % wavenumber f = -Nyq_f : df : Nyq_f-df; % frequency data_wavenumberfrequency = zeros(size(data_spacetime)); % initialize data % Compute 2D Discrete Fourier Transform for i1 = 1 : Nx for j1 = 1 : Nt for i2 = 1 : Nx for j2 = 1 : Nt data_wavenumberfrequency(j1,i1) = data_wavenumberfrequency(j1,i1) + ... data_spacetime(j2,i2)*exp(-1i*(2*pi)*(k(i1)*x(i2)+f(j1)*t(j2)))*dx*dt; end end end end figure; subplot(3,1,1); imagesc(k,f,abs(data_wavenumberfrequency)); colorbar; v = caxis; title('2D DFT'); fft2result = fftshift(fft2(data_spacetime))*dx*dt; subplot(3,1,2); imagesc(k,f,abs(fft2result)); colorbar; caxis(v); title('FFT2'); subplot(3,1,3); imagesc(k,f,abs(fft2result-data_wavenumberfrequency)); colorbar; caxis(v); title('Difference');

【本文地址】

公司简介

联系我们