function varargout = xdevdemo304(WformIn1,number)
%XDEVDEMO304 Recovers the primary sinewaves of the waveform.
%
% WFORMOUT = XDEVDEMO304(WFORMIN1,NUMBEROFTONES)finds the local peaks
% in the fft of the input data WFORMIN1 assuming the highest peak is 
% the first peak. The rebuilt waveform, using only those frequencies, is 
% then returned.  NUMBEROFTONES is the maximum number of tones to use in
% reconstructing the waveform.  If the maximum number of tones found is 
% less than NUMBEROFTONES then only the tones found are used.
%
% Example:
%   WformOut = xdevdemo304(WformIn1,10);
%
% See also DEMOS, XDEVDEMO, XDEVDEMO300.
%

%   $Author: Ldean $Revision: 1.0 $  $Date: 9/24/02 5:32p $

%   Testing:
%   This function was tested on the WAVEMASTER series scope.
%

% $Notes: 
% $EndNotes 

% $Description:
% $
if nargin==0
    help(mfilename)
    return;
end;

if nargin==1
	number = 1;
end;

% Don't let people choose more than 20 harmonics
n=min(20,number);

% Prealocate memory.
WformOut=zeros(size(WformIn1));
for c=1:n
	x=localxdevdemo301(WformIn1); % Get primary wave from current data set
	WformIn1=WformIn1-x;
	WformOut=WformOut+x;
end;
varargout{1} = WformOut;

function WformOut = localxdevdemo301(WformIn1)
% Take the FFT of the data and find the largest signal above DC
% then generate a sin wave with this frequency.

% Get the length of the signal
N=length(WformIn1);

% Take the FFT of the data removing DC
fftdata = fft(WformIn1-mean(WformIn1));

% Find the maximum peak in the FFT
[y,idxTone] = max(fftdata);
idxTone=idxTone;
% Calculate the frequency that this peak represents.
freq=(idxTone-1)/N;

% Generate the sin wave that represents this frequency using the angle from
% the FFT data.
a1=angle(fftdata(idxTone));
a2=angle(fftdata(N-idxTone+2));
if a1>=0
    if a2>0
        disp('Imaginary Results in FFT.  We have a problem.');
        ang = 0;        
    else
        ang = a1+pi/2; % cos(ang)
    end;
else
    if a2>=0
        ang = a1+pi/2; %sin(ang)
    else
        disp('Imaginary Results in FFT.  We have a problem.');
        ang = 0;        
    end;
end;
WformOut=(2*abs(y)/N)*sin(2*pi*[0:N-1]'*freq+ang);

