✅作者简介：热爱科研的Matlab仿真开发者，修心和技术同步精进，matlab项目合作可私信。
🍎个人主页：Matlab科研工作室
🍊个人信条：格物致知。
更多Matlab仿真内容点击👇
⛄ 内容介绍
快速傅里叶变换(Fast Fourier Tranformation,FFT)是将一个大点数N的DFT分解为若干小点的DFT的组合,将用运算工作量明显降低,从而大大提高了离散傅里叶变换(DFT)的计算速度.因各个科学技术领域广泛的使用了FFT技术,它大大推动了信号处理技术的进步,现已成为数字信号处理强有力的工具.本论文将比较全面地叙述各种快速傅里叶变换算法原理,特点,并完成了基于MATLAB的实现.
⛄ 部分代码
% Demoedft.m demonstrates performance of function EDFT by iterations
% and provide comparison with Nonuniform DFT (Matlab NUFFT) output.
% EDFT call line is [F,S,Stopit]=edft(X,N,1,W,tk),
% where DEMO data X are calculated as sum of the following components:
% - real sinusoid at normalized frequency 0.137, magnitude 2,
% - complex exponent at frequency 0.137+1/(4*K), magnitude 1, where K=64
% is the length of sequence X,
% - mean value (frequency 0), magnitude 1,
% - pulse in frequency ranges [-0.137 0],
% - noise generated by rand function (SNR~35 dB).
% Number of frequencies N=512 (the length of DFT).
% Number of samples in sequence with missing data KK=32.
% So, the true spectrum of DEMO sequence (without noise part) is ploted
% in red color and looks like greeting 'HI', where 'I' actually are two
% close peaks at normalized frequencies 0.137 and 0.137+1/(4*K).
% In the first, demoedft.m plot real parts of the input sequences:
% - Subplot(511) display DEMO data X1 consisting of 64 uniform samples;
% - Subplot(512) outputs an example of missing data X2, where 32 randomly
% taken samples in X1 replaced by NaN;
% - Subplot(513) shows gapped data X3, where two gaps by 16 NaN each
% inserted in X1.
% - Subplot(514) outputs nonuniformly sampled DEMO data X4,
% tkn2=[1:K]+0.5*(rand(1,K)-0.5);
% - Subplot(515) showing sparse (32 samples) nonuniformly sampled data X5,
% tkn3=tkn2(1:2:K);
% NaN are not displayed in sublopts 512 and 513.
% To proceed you should select one of these sequences displayed in figure.
% In the next, fifteen (I=15) NEDFT iterations are performed. You will be
% asked to strike any key after 1,2,3,5 and 10 iteration. Plots in blue
% colour are equal to NUFFT and can be used for comparison with EDFT
% output (green). To resolve two close sinusoids (at 1/4 of Nyquist), EDFT
% should increase resolution at least 4 times in comparison with NUFFT.
%
% EDFT and NUFFT outputs are displayed in four subplots:
% -Subplot(221) show Power Spectral Density as 10*log10(abs(F).^2/N)) in
% normalized frequencies range [-0.5 0.5[.
% -Subplot(222) display Power Spectrum as 10*log10(abs(S).^2). The plot
% indicate power of sinusoids in X.
% -Subplot(413) plot division F./S/K titled as Relative frequency resolution
% and expose how the frequency resolution of EDFT algorithm is changed
% along the frequency axes in respect to NUFFT. The first iteration
% showing that NUFFT and EDFT have almost equal resolution on all
% frequencies in range [-0.5 0.5[. During the next EDFT iterations we
% can discover that the relationship F./S/K for EDFT analysis is not that
% simple. Although sum(F./S/K)=N and remains constant for all iterations,
% EDFT have ability to increase resolution around the powerful narrowband
% components (sinusoids) and decrease resolution at frequency range where
% data X have weak power components. EDFT is called as high resolution
% method and that's true, but with the following remark- EDFT keeping
% the same 'summary' resolution as NUFFT has or, in other words, squares
% under curves for NUFFT (blue) and EDFT (green) are equal. The maximum
% frequency resolution is limited by value of division N/K. For example,
% if K=64 and N=512, then EDFT can potentially improve frequency
% resolution 512/64=8 times.
% -Subplot(414) plot real parts of reconstructed sequences obtained in
% result of Inverse Fast Fourier Transform to outputs of NUFFT and EDFT,
% correspondingly. It is well known that reconstruction of data by
% applying ifft(fft(X,N)) for N>length(X), return sequence where initial
% data X are padded with zeros to length N. That's true also for NUFFT
% and for the first iteration of EDFT. The next iterations are showing
% an ability of EDFT to obtain reconstructed data consisting of input
% sequence X values plus non-zero extrapolation of X to length N. The
% sequence X is extrapoated in both directions - backward and forward.
% NaNs in sequences [X2] and [X3] are interpolated and replaced with
% reconstructed values, while for nonuniform sequences [X4] and [X5],
% the reconstructed data are re-sampled at uniform time points tn=1:N.
%
% References:
% 1) Vilnis Liepins, A spectral estimation method of nonuniformly sampled
% band-limited signals.Automatic Control and Computer Sciences,Vol.28,
% No.2, pp.66-77, 1994.
% 2) Vilnis Liepins, An algorithm for evaluation a discrete Fourier
% transform for incomplete data, Automatic control and computer sciences,
% Vol.30, No.3, pp.27-40, 1996.
clear
disp('EDFT DEMO program started...');
% Set the length of input data (K), DFT (N), number of iterations (I),
% relative central frequency (fc) in range ]0,0.5] and SNR:
K=64;
N=512;
I=15;
fc=0.137;
SNR=35;
% Uniformly and nonuniformly spaced time and frequency vectors
tk=0:K-1; % uniform time vector (K samples, sampling period T=1)
tn=0:N-1; % uniform time vector (N samples, sampling period T=1)
tkn2=(1:K)+0.5*(rand(1,K)-0.5); % nonuniform (mean sampling period T=1)
tkn3=tkn2(1:2:K); % sparse nonuniform (mean sampling period T=2)
fn=(-ceil((N-1)/2):floor((N-1)/2))/N; % normalized frequencies [0.5,0.5[
ti=tk-K/2-0.5; % uniform centerred time vector for pulse
tin2=tkn2-K/2-0.5; % nonuniform centerred time for pulse
% Generate X1 - uniform sequence (K samples)
Xph=2*pi*rand(1,2);
XC1=2*sin(2*pi*fc*tk+Xph(1,1)); % sinusoid
XC2=exp(1i*(2*pi*(fc+1/4/K)*tk+Xph(1,2))); % complex exponent
XC3=sin(pi*fc*ti)./(ti+eps).*exp(-1i*pi*fc*ti); % rectangular pulse
XC4=ones(1,K); % mean value
X1=XC1+XC2+XC3+XC4; % signal w/o noise
sigma=sqrt((2+1+1)/10^(SNR/10)); % sigma for given SNR
XC5=sigma*(randn(1,K)+1i*randn(1,K))/sqrt(2);% complex noise
SNR1=round(10*log10((X1*X1')/(XC5*XC5')));% calculated SNR
X1=X1+XC5; % X1 signal+noise (EDFT and NUFFT input)
% Generate X2 - nonuniform sequence with NaN: (K/2 samples, K/2 NaN)
MNaN=floor(K/2); % MNaN- number of NaN in X
KK=K-MNaN; % KK- length of X2,X3 without NaN
X2=X1;XX2=[];t2=[];
KNaN=(1:2:MNaN*2)+round(rand(1,MNaN)); % KNaN indicate NaN in sequence
X2(KNaN)=NaN*ones(1,MNaN); % sequence with NaN (EDFT input)
XX2(1:KK)=X1(~isnan(X2)); % sequence (NUFFT input)
t2(1:KK)=tk(~isnan(X2)); % time vector (NUFFT input)
% Generate X3 - gapped sequence (K/2 samples, 2 x K/4 NaN)
X3=X1;XX3=[];t3=[];k4=floor(K/4);
ki=round([rand*k4+(1:k4) MNaN+rand*(MNaN-k4)+(1:(MNaN-k4))]);
X3(ki)=NaN*ones(1,MNaN); % gapped data with NaN (EDFT input)
XX3(1:KK)=X1(~isnan(X3)); % sequence (NUFFT input)
t3(1:KK)=tk(~isnan(X3)); % time vector (NUFFT input)
% Generate X4 - nonuniform sequence (K samples)
XN1=2*sin(2*pi*fc*tkn2+Xph(1,1)); % sinusoid
XN2=exp(1i*(2*pi*(fc+1/4/K)*tkn2+Xph(1,2)));% complex exponent
XN3=ones(1,K); % mean value
XN4=sin(pi*fc*tin2)./(tin2+eps).*exp(-1i*pi*fc*tin2); % rectangular pulse
X4=XN1+XN2+XN3+XN4; % signal w/o noise
SNR2=10*log10((X4*X4')/(XC5*XC5'));% calculated SNR
X4=X4+XC5; % signal+noise
% Generate X5 - sparse nonuniform sequence (half of X4 samples)
X5=X4(1:2:K);K5=length(X5);
% Plot real part of sequences X1, X2, X3, X4 and X5.
figure(1)
clf
Xmin=min(real([X1 X4]))-1;
Xmax=max(real([X1 X4]))+1;
subplot(511)
stem(tk+1,real(X1))
axis([0 K+1 Xmin Xmax])
xlabel('Sample number')
ylabel('real(X1)')
title(['Input sequence X1 - uniform data: ',int2str(K),' samples, SNR=' ...
,int2str(round(SNR1))])
subplot(512)
stem(tk+1,real(X2))
axis([0 K+1 Xmin Xmax])
xlabel('Sample number')
ylabel('real(X2)')
title(['Input sequence X2 - missing data: ',int2str(KK),' samples and ' ...
,int2str(MNaN),' NaN'])
subplot(513)
stem(tk+1,real(X3))
axis([0 K+1 Xmin Xmax])
xlabel('Sample number')
ylabel('real(X3)')
title(['Input sequence X3 - gapped data: ',int2str(KK),' samples and ' ...
,int2str(k4),'+',int2str(MNaN-k4),' NaN'])
subplot(514)
stem(tkn2,real(X4))
axis([0 K+1 Xmin Xmax])
xlabel('Sample number')
ylabel('real(X4)')
title(['Input sequence X4 - nonuniform data: ',int2str(K), ...
' samples, SNR=',int2str(round(SNR2))])
subplot(515)
stem(tkn3,real(X5))
axis([0 K+1 Xmin Xmax])
xlabel('Sample number')
ylabel('real(X5)')
title(['Input sequence X5 - sparse nonuniform data: ' ...
,int2str(K/2),' samples'])
% Selecting of sequence for NEDFT input and calculating DFT output
disp('DEMO data are calculated as sum of the following components:');
disp([' - real sinusoid at normalized frequency ',num2str(fc), ...
', magnitude 2;']);
disp([' - complex exponent at frequency ',num2str(fc),'+1/(4*' ...
,int2str(K),'), magnitude 1;']);
disp(' - mean value (frequency 0), magnitude 1;');
disp([' - rectangular pulse in frequency ranges [-',num2str(fc),' 0];']);
disp([' - noise generated by randn function (SNR~',int2str(SNR),' dB).']);
DEMO=1;
while DEMO==1||DEMO==2||DEMO==3||DEMO==4||DEMO==5
disp(' ');
disp('Input [1] to select uniform sequence X1');
disp('Input [2] to select missing sample sequence X2');
disp('Input [3] to select gapped sequence X3');
disp('Input [4] to select nonuniformly sampled sequence X4');
disp('Input [5] to select sparse nonuniformly sampled sequence X5');
DEMO=input('Input sequence # or any other valid key to end of DEMO:');
if DEMO==1
X=X1;XX=X1;tkk=tk;KP=K;h=2; % X1 - uniform data
elseif DEMO==2
X=X2;XX=XX2;tkk=t2;KP=KK;h=3; % X2 - missing data
elseif DEMO==3
X=X3;XX=XX3;tkk=t3;KP=KK;h=4; % X3 - gapped data
elseif DEMO==4
X=X4;XX=X4;tkk=tkn2;KP=K;h=5; % X4 - nonuniform data (K samples)
elseif DEMO==5
X=X5;XX=X5;tkk=tkn3;KP=K5;h=6; % X5 - sparse nonuniform data
else
disp('...end of EDFT DEMO.') % End DEMO if other valid key entered
return
end
% Calculate outputs for NUFFT used in 4 plottings.
dft_out=nufft(XX,tkk,fn);
sub1_PWD=10*log10(abs(dft_out).^2/N);
sub2_PS=20*log10(abs(dft_out/KP));
sub3_FR=ones(1,N)*KP/K;
sub4_RD=real(ifft(ifftshift(dft_out)));
% Set values for input argument W used for the first EDFT iteration.
W=ones(1,N);
% Get EDFT outputs F and S for iteration it.
for it=1:I
if DEMO==4||DEMO==5
[F,S,Stopit]=edft(X,N,1,W,tkk); %Basic algorithm
else
[F,S,Stopit]=edft(X,N,1,W); %Faster algorithm
F=fftshift(F); %Adjust size
S=fftshift(S); %of output
end
% Calculate outputs for EDFT iteration it used in 4 plottings.
sub1=10*log10(abs(F).^2/N);
sub2=20*log10(abs(S));
sub3=real(F./S/K);
sub4=real(ifft(ifftshift(F)));
frsig=[-fc 0 fc fc+1/4/K];
% Plots Power Spectral Densities in subplot221.
figure(h)
clf
subplot(221)
plot([frsig;frsig],[10*log10(N)*[1 1 1 1];-100*[1 1 1 1]],'r:');
hold on
plot([-fc;0],[10*log10(pi^2/N);10*log10(pi^2/N)],'r:');
plot(fn,sub1,'g-',fn,sub1_PWD,'b-')
hold off
axis([-0.5 0.5 -100 50])
xlabel('Normalized frequency')
ylabel('10*log[abs(F).^2/length(F)]')
title(['Power Spectral Density (iter.',int2str(it),')'])
% Plots Power Spectrums in subplot222.
subplot(222)
plot([frsig;frsig],[[0 0 0 0];-100*[1 1 1 1]],'r:');
hold on
plot([-fc;0],[20*log10(pi/K);20*log10(pi/K)],'r:');
plot(fn,sub2,'g-',fn,sub2_PS,'b-')
hold off
axis([-0.5 0.5 -100 50])
xlabel('Normalized frequency')
ylabel('10*log[abs(S).^2]')
title(['Power Spectrum (iter.',int2str(it),')'])
% Plots Frequency Resolutions in subplot413.
subplot(413)
plot([frsig;frsig],[[0 0 0 0];N/K*[1 1 1 1]],'r:');
hold on
plot(fn,sub3,'g-',fn,sub3_FR,'b-')
plot([-fc;0],[1;1],'r:');
hold off
axis([-0.5 0.5 0 N/K])
xlabel('Normalized frequency')
ylabel(['(F./S)/length(X',int2str(DEMO),')'])
title(['Relative frequency resolution (iter.',int2str(it),')'])
% Plots orriginal (red) and Reconstructed Data in subplot414.
subplot(414)
plot(tk,real(X1),'r-',tn,sub4,'g-',tn,sub4_RD,'b-')
axis([1 N Xmin Xmax])
xlabel('Sample number')
ylabel('real(ifft(F))')
title(['Sequence X',int2str(DEMO),': True (red) and reconstructed' ...
' by EDFT (green, iter.',int2str(it),') and NUFFT (blue)'])
% Waiting for keyboard action
drawnow
if it==1
disp(' ');
disp(['Plots True Spectrum [red color], EDFT [green colour] and' ...
' NUFFT [blue colour].']);
end
if it==1||it==2||it==3||it==5||it==10
disp(['Calculate [F,S]=edft(X',int2str(DEMO),',',int2str(N),',' ...
'',int2str(it),',...) and ifft(F). Strike any key to continue.']);
pause
end
% Set EDFT input argument W for next iteration.
if DEMO==4||DEMO==5
W=S;
else
W=ifftshift(S); %Adjust size of output
end
end
disp(['Calculate [F,S]=edft(X',int2str(DEMO),',' ...
'',int2str(N),',',int2str(it),',...) and ifft(F).']);
end
⛄ 运行结果
⛄ 参考文献
[1]谭子尤, 张雅彬. 离散傅里叶变换快速算法的研究与MATLAB算法实现[J]. 中国科技信息, 2006(22):3.
[2]马维祯. 多维离散傅里叶变换DFT(2n;k)算法[C]// 中国第十届电路与系统学术年会. 1992.
