demoedft.m

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

% Demoedft.m demonstrates performance of function EDFT by iterations and provide comparison with 
% traditional discrete Fourier transform (DFT) output. EDFT call line is [F,S,Stopit]=edft(X,N,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 (length of Fourier transform) N=1024.
% Number of samples (length of input sequence) 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, demoedft.m plot real parts of the input sequences:
% - Subplot(311) display DEMO data X1 consisting of 64 samples;
% - Subplot(312) outputs an example of sparse data X2, where 32 randomly taken samples in X1 replaced by NaN;
% - Subplot(313) shows gapped data X3, where two gaps by 16 NaN each inserted in X1.
% NaN are not displayed in sublopts 312 and 313. 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 only 32 available samples
% will be processed by DFT and EDFT algorithms. 
%
% In the next, fifteen (I=15) EDFT iterations are performed starting with edft(X,N,1,ones(1,N))=fft(X,N).
% You will be asked to strike any key after 1,2,3,5 and 10 iteration. Plots in blue colour are equal
% to traditional DFT analysis and can be used for comparison with EDFT outputs (green colour). To resolve
% two close sinusoids (at 1/5 of Nyquist), EDFT should increase resolution at least 5 times 
% (in comparison with traditional DFT, see Subplot(413)).
%
% EDFT 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[.
% Note, that results of the first iteration coincide to periodogram: 10*lg(abs(fft(X,N).^2/N)).
% - Subplot(222) display Power Spectrum as 20*log10(abs(S)). That plot indicates the power of sinusoids in data 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 DFT analysis. The first iteration plot
% F./S/K=const and showing that traditional DFT have equal resolution for all frequencies in range [-0.5 0.5[.
% Next EDFT iterations discover that the relationship F./S/K for EDFT analysis is not that simple as for
% traditional DFT. Although sum(F./S/K)=N and remains constant for all iterations, EDFT have ability to
% increase resolution (in plot appears values >1) around the powerful narrowband components (sinusoids) 
% and decrease resolution (in plot appears values <1) at frequencies 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 traditional DFT or, in other words, squares under curves (F./S)/K for traditional
% DFT (blue colour) and EDFT (green colour) analysis are equal. Maximum frequency resolution is limited by
% value of division N/K. For example, if K=64 and N=1024 then EDFT 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 EDFT - ifft(dft_out) and ifft(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 padded with zeros to 
% length N. That result is displayed for the first iteration. The next EDFT iterations provide ability to obtain
% reconstructed sequence consisting of initial data X values plus non-zero extrapolation of X to length N. As shown
% in subplot 414 data X are extrapolated in both directions - backward and forward. But in case of NaN, data also
% are interpolated and 32 NaN are replaced with reconstructed values.
%
% See attached reference article for additional info. 	
%
% Last modified by Vilnis Liepins on 1/3/2008.
% E-mail:	[e-mail]
% Reference: V. Ya. Liepin'sh, "An algorithm for evaluation a discrete Fourier transform for incomplete data",
% Automatic control and computer sciences, Vol.30, No.3, pp.27-40, 1996 (ISSN0132-4160).

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

% Set length for input data 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);

% uniformly spaced time and frequency vectors
	tk=1:K;				% uniform time vector (K samples, sampling period T=1)
	tn=1:N;				% uniform time vector (N samples, sampling period T=1)
	fn=-0.5:1/N:0.5-1/N;		% uniform normalized frequency set (sampling frequency- 1)

% Generate X1 - uniformly sampled EDFT input sequence (K samples)
	% signal
	X1=2*sin(2*pi*fc*tk)+exp(i*2*pi*(fc+1/5/K)*tk)...
	+sin(pi*fc*(tk-K/2-0.5))./(tk-K/2-0.5).*exp(-i*pi*fc*tk)+1;
	% add noise
	X1=X1+0.1*(rand(size(tk))-0.5)+i*0.1*(rand(size(tk))-0.5);
	
% Generate X2 - nonuniformly sampled EDFT input: sparse sequence (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;k=[];
	while length(k)~=MNaN,
	    k=find(sprandn(K,1,0.8));	% k- indicate NaN in input sequence 
	end
        X2(k)=zeros(1,MNaN);XX2=X2;	% XX2- X2 where NaN replaced with zeros (FFT input) 
	X2(k)=NaN*ones(1,MNaN);		% X2 - sparse data with NaN (EDFT input)
	
% Generate X3 - nonuniformly sampled EDFT input: gapped sequence (K/2 samples, 2 x K/4 NaN)

	X3=X1;
	k=round([rand*K/4+(1:K/4) K/2+rand*K/4+(1:K/4)]);  
        X3(k)=zeros(1,MNaN);XX3=X3;	% XX3- X3 where NaN replaced with zeros (FFT input)
	X3(k)=NaN*ones(1,MNaN);		% X3 - gapped data with NaN (EDFT input)
	

% Plot real part of sequences X1, X2 and X3.

	figure
	Xmin=min(real(X1))-1;
	Xmax=max(real(X1))+1;
	subplot(311)
	stem(tk,real(X1))
	axis([0 K+1 Xmin Xmax])
	xlabel('Sample number')
	ylabel('real(X1)')
	title(['Input sequence X1 - uniform data: ',int2str(K),' samples'])
	subplot(312)
	stem(tk,real(X2))
	axis([0 K+1 Xmin Xmax])
	xlabel('Sample number')
	ylabel('real(X2)')
	title(['Input sequence X2 - Sparse data: ',int2str(KK),' samples and ',int2str(MNaN),' NaN'])
	subplot(313)
	stem(tk,real(X3))
	axis([0 K+1 Xmin Xmax])
	xlabel('Sample number')
	ylabel('real(X3)')
	title(['Input sequence X3 - Gapped data: ',int2str(KK),' samples and 2x',int2str(K/4),' NaN'])

% Selecting of sequence for EDFT input. Input DEMO==2 or 3 choose sequences X2 or X3 where
% K/2 of K samples replaced with NaN to demonstrate EDFT performance for incomplete data.

	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;XX=X1;KK=K;end	% X=X1 - uniform data
        if DEMO==2,X=X2;XX=XX2;end		% X=X2 - sparse data		
	if DEMO==3,X=X3;XX=XX3;end		% X=X3 - gapped data
	else
	disp('...end of DEMO')			% End DEMO if other symbol entered
	return
	end

% Calculate outputs for DFT (periodogram) used in 4 plottings.
	
	dft_out=fft(XX,N);
	sub1_PWD=10*log10(fftshift(abs(dft_out).^2/N));
	sub2_PS=10*log10(fftshift(abs(dft_out/KK).^2));
	sub3_FR=ones(1,N)*KK/K;
	sub4_RD=real(ifft(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,

	[F,S,Stopit]=edft(X,N,1,W);

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

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

	sub1=10*log10(fftshift(abs(F).^2/N));
	sub2=10*log10(fftshift(abs(S).^2));
	sub3=real(fftshift(F./S)/K);
	sub4=real(ifft(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

% Set EDFT input argument W for next iteration.

	W=S;

	end

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