demonedft.m

程序之所以称之为扩展功能的DFT

% Demonedft.m demonstrates performance of NEDFT algorithm by iterations in case of nonuniform sampling and
% provide comparison with traditional nonuniform discrete Fourier transform (DFT) output. NEDFT call line
% is [F,S,Stopit]=nedft(X,tk,fn,1,W), 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/(5*K), magnitude 1 (K=64),
% - rectangular pulse in frequency ranges [-0.137 0],
% - mean value (frequency 0), magnitude 1 and
% - noise generated by rand function (SNR=30 dB).
% Number of frequencies in fn (length of Fourier transform) N=1024.
% Number of samples (length of input sequence X) K=64 or 32.
% So, the true spectrum of DEMO sequence (without noise part) looks like greeting 'HI', where 'I' actually
% are two close peaks at normalized frequencies 0.137 and 0.137+1/(5*K) - 1/5 of Nyquist.
%
% In the first, demonedft.m plots real parts of the input sequences:
% - Subplot(311) display DEMO data X1 consisting of 64 samples, tk=[1:K];
% - Subplot(312) outputs nonuniformly sampled DEMO data X2, tkn2=[1:K]+0.5*(rand(1,K)-0.5);
% - Subplot(313) showing sparse (32 samples) nonuniformly sampled data X3, tkn3=tkn2(1:2:K); 
% To proceed you should select one of these sequences - enter [1] to see DEMO for 64 uniformly sampled
% data or select [2] or [3], where 64 or 32 nonuniformly sampled data will be processed by DFT and NEDFT
% algorithms.	
%
% 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 DFT analysis and can be used 
% for comparison with NEDFT outputs (green colour). To resolve two close sinusoids (at 1/5 of Nyquist),
% NEDFT should increase resolution at least 5 times (in comparison with DFT, see Subplot(413)).
%
% NEDFT and DFT 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 20*log10(abs(S)). 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 NEDFT algorithm is changed along the frequency axes in respect to DFT analysis. The first iteration showing
% that DFT and NEDFT analysis have almost equal resolution for all frequencies in range [-0.5 0.5[. During the next
% NEDFT iterations we can discover that the relationship F./S/K for NEDFT analysis is not that simple as for DFT.
% Although sum(F./S/K)=N and remains constant for all iterations, NEDFT have ability to increase resolution 
% around the powerful narrowband components (sinusoids) and decrease resolution at frequencies where data X have
% weak power components. NEDFT is called as high resolution method and that's true, but with the following remark-
% NEDFT keeping the same 'summary' resolution as DFT or, in other words, squares under curves for DFT (blue colour)
% and NEDFT (green colour) analysis are equal. The maximum frequency resolution is limited by value of division N/K.
% For example, if K=64 and N=1024, then NEDFT can potentially improve frequency resolution 1024/64=16 times. 
% - Subplot(414) plot real parts of reconstructed sequences obtained as result of inverse fast Fourier transform
% to outputs of DFT and NEDFT - ifft(fftshift(dft_out)) and ifft(fftshift(F)), correspondingly. It is well known
% that reconstruction of data by applying ifft(fft(X,N)) for length(X)<N, return sequence where initial data X are
% padded with zeros to length N. That result is true also for DFT and for the first iteration of NEDFT. The next
% iterations are showing ability of NEDFT to obtain reconstructed data consisting of input sequence X values plus
% non-zero extrapolation of X to length N. As shown in subplot 414, sequence X are extrapolated in both directions
% - backward and forward. Also for nonuniformly sampled input sequences [X2] and [X3], the reconstructed data are
% re-sampled at uniform time set tn=1:N.
%
% Last modified by Vilnis Liepins on 1/3/2008.
% E-mail:	[e-mail]

	clear
	disp('NEDFT DEMO program started...');

% Set the length for input sequence X (K), DFT (N) and number of iterations (I).
	
	K=64;
	N=1024;
	I=15;

% Relative central frequency of the test signal (value in range ]0,0.5]).
	fc=0.137;
			
% Generate always the same noice for input sequences.
	%rand('seed',0);

% Generate uniformly and nonuniformly spaced time and frequency vectors.

	tk1=1:K;			% uniform time vector (K samples, sampling period T=1)
	tkn2=[1:K]+0.5*(rand(1,K)-0.5);	% nonuniform time vector (K samples, mean sampling period Ts=1)
	tkn3=tkn2(1:2:K); 		% sparse nonuniform time vector (K/2 samples, mean sampling period Ts=2)
	fn=-0.5:1/N:0.5-1/N;		% uniform normalized frequency set (sampling frequency- 1)
	tn=1:N;				% uniform time vector (N samples, sampling period T=1)

% Generate uniformly sample input sequence X1 (K samples).
	% signal
	X1=2*sin(2*pi*fc*tk1)+exp(i*2*pi*(fc+1/5/K)*tk1)...
	+sin(pi*fc*(tk1-K/2-0.5))./(tk1-K/2-0.5).*exp(-i*pi*fc*tk1)+1;
	% noice
	X1=X1+0.1*(rand(size(tk1))-0.5)+i*0.1*(rand(size(tk1))-0.5);


% Generate nonuniformly sampled sequence X2 (K samples)
        % signal
	X2=2*sin(2*pi*fc*tkn2)+exp(i*2*pi*(fc+1/5/K)*tkn2)...
	+sin(pi*fc*(tkn2-K/2))./(tkn2-K/2).*exp(-i*pi*fc*tkn2)+1;
	% noice
	X2=X2+0.1*(rand(size(tkn2))-0.5)+i*0.1*(rand(size(tkn2))-0.5);

% Generate nonuniformly sampled sequence X3 (taking half of X2 samples)
	
	X3=X2(1:2:K);

% Plot the real part of data X1, X2 and X3.

	figure
	Xmin=min(real(X1))-1;
	Xmax=max(real(X1))+1;
	subplot(311)
	stem(tk1,real(X1))
	axis([0 K+1 Xmin Xmax])
	xlabel('Sample number')
	ylabel('real(X1)')
	title(['Input sequence X1 - uniform sampling: ',int2str(K),' samples'])
	subplot(312)
	stem(tkn2,real(X2))
	axis([0 K+1 Xmin Xmax])
        xlabel('Sample number')
	ylabel('real(X2)')
	title(['Input sequence X2 - nonuniform samling: ',int2str(K),' samples'])
	subplot(313)
	stem(tkn3,real(X3))
	axis([0 K+1 Xmin Xmax])
        xlabel('Sample number')
	ylabel('real(X3)')
	title(['Input sequence X3 - sparse nonuniform sampling: ',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/(5*' int2str(K) '), magnitude 1']);
	disp(['   - rectangular pulse in frequency ranges [-' num2str(fc) ' 0],']);
        disp('   - mean value (frequency 0), magnitude 1 and');
        disp('   - noise generated by rand function (SNR~30 dB).');
	disp('Press Enter or input [1] to select uniformly sampled sequence X1');
	disp('Input [2] to select sparse sequence X2');
	disp('Input [3] to select gapped sequence X3');
	DEMO=input('Select input sequence for DEMO:');  
	if isempty(DEMO)|DEMO==1|DEMO==2|DEMO==3,
	if isempty(DEMO)|DEMO==1,
	    X=X1;tk=tk1;KK=K;			% X=X1 - uniform data (K samples)
	    E=exp(-i*2*pi*tk1.'*fn);		% E - complex exponents matrix
	    dft_out=X*E; 			% Calculate DFT output
	end		
        if DEMO==2,
	    X=X2;tk=tkn2;KK=K;			% X=X2 - nonuniform data (K samples)
	    E=exp(-i*2*pi*tkn2.'*fn);		% E - complex exponents matrix
	    dft_out=X*E; 			% Calculate DFT output
	end							
	if DEMO==3,
	    X=X3;tk=tkn3;KK=K/2;		% X=X3 - sparse nonuniform data (K/2 samples)
	    E=exp(-i*2*pi*tkn3.'*fn);		% E - complex exponents matrix
	    dft_out=X*E; 			% Calculate DFT output
	end		
	else
	   disp('...end of DEMO')		% End DEMO if other symbol entered
	   return
	end

% Calculate outputs for DFT used in 4 plottings.

	sub1_PWD=10*log10(abs(dft_out).^2/N);
	sub2_PS=10*log10(abs(dft_out/KK).^2);
	sub3_FR=ones(1,N)*KK/K;
	sub4_RD=real(ifft(fftshift(dft_out)));


% Set values for input argument W used for the first NEDFT iteration.

	W=ones(1,N);

% Get NEDFT outputs F and S for iteration it.

	for it=1:I,

	[F,S,Stopit]=nedft(X,tk,fn,1,W);

% Break NEDFT iterations if get ill conditioned matrix or results may be inaccurate.
	if Stopit(2)==1|Stopit(2)==2,break,end 


% Calculate outputs for NEDFT iteration it used in 4 plottings.

	sub1=10*log10(abs(F).^2/N);
	sub2=10*log10(abs(S).^2);
	sub3=real(F./S/K);
	sub4=real(ifft(fftshift(F)));

% Plots Power Spectral Densities in subplot221.
	
	clf
	subplot(221)
	plot(fn,sub1,'g-',fn,sub1_PWD,'b-')
	axis([-0.5 0.5 -100 50])
	xlabel('Normalized frequency')
	ylabel('10*log[abs(F).^2/length(F)]')
	title(['Power Spectral Density (it.',int2str(it),')'])
	
% Plots Power Spectrums in subplot222.

	subplot(222)
        plot(fn,sub2,'g-',fn,sub2_PS,'b-')
	axis([-0.5 0.5 -100 50])
	xlabel('Normalized frequency')
	ylabel('10*log[abs(S).^2]')
	title(['Power Spectrum (it.',int2str(it),')'])

% Plots Frequency Resolutions in subplot413.

	subplot(413)
	plot(fn,sub3,'g-',fn,sub3_FR,'b-')
	axis([-0.5 0.5 0 N/K])
	xlabel('Normalized frequency')
       	ylabel('(F./S)/length(X)')
	title(['Relative frequency resolution (it.',int2str(it),')'])

% Plots Reconstructed data in subplot414.

	subplot(414)
	plot(tn,sub4,'g-',tn,sub4_RD,'b-')
	axis([1 N Xmin Xmax])
	xlabel('Sample number')
	ylabel('real(ifft(F))')
	title(['Reconstructed sequence (it.',int2str(it),')'])

% Waiting for keyboard action

	drawnow
        if it==1|it==2|it==3|it==5|it==10,
	    disp('Strike any key to continue.')
	    pause
	end

% Calculate NEDFT input argument W for next iteration.

	W=S;

	end

	disp('...end of DEMO');