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

参考帖子
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');