function varargout = xdevdemo303(WformIn1,number)
%XDEVDEMO303 Recovers the primary sinewaves of the waveform.
%
%  WFORMOUT = XDEVDEMO303(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 = xdevdemo303(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;


N=length(WformIn1);
meanValue = mean(WformIn1);
[val,freqidx]=localpeaksInternal(WformIn1-meanValue); % Shown below

if isempty(val)
	WformOut = zeros(size(WformIn1));
else
    numTones = length(val);
    if nargin==2
        if (numTones>number)
            useTone=[1:number,numTones-number+1:numTones];
            val=val(useTone);
            freqidx=freqidx(useTone);
        end;
    end;
    freqidx = freqidx-1;
    idxImage = find(freqidx>floor(N/2));
    val(idxImage)=-1*val(idxImage);

    freq=(2*pi*[0:N-1]'/N);
    WformOut=zeros(size(WformIn1));
    
    idx=min(length(val),20); % This is just to make it faster (DEMO Code).
    
    for c=1:idx;
        WformOut=WformOut+val(c)*sin(freq*freqidx(c));
    end;
end;	
varargout{1} = WformOut;

function [localPeakValue, localPeakIdx] = localpeaksInternal(waveformin)
% LOCALPEAKS	Find the amplitude of the primary frequecies.
%	
% This function will find the local peaks in the fft of the data using
% the given window size.  It does not report the peak of the DC signal 
% and assumes the highest peak is the first peak.
%
%
% Notes:
% End of Notes.

% Calculate FFT and frequency spacing 
N = length(waveformin);
fftdata = abs(fft(waveformin))/N;

% We are assuming that the first peak is the maximum except for DC
%	we are removing the DC component.

% If a window width is not specified we will use half the
% distance from DC to first Peak as the window width.
[maxPeak,idxPeak] = max(fftdata);
windowWidth = floor(idxPeak/2)+1;

% Find all the peak frequencies in the given window
M = floor(N/windowWidth);
fdata = reshape(fftdata(1:M*windowWidth),windowWidth,M);
[localPeakValue,idxLocalPeak]=max(fdata);
idxLocalPeak = idxLocalPeak(2:end)+([1:M-1]*windowWidth);
localPeakValue = localPeakValue(2:end);

% Find the noise floor
noisefloor = mean(fftdata);

% Only report back frequencies with values above 2*noise floor
[newLocalPeakidx] = find(localPeakValue>2*noisefloor);
localPeakIdx = idxLocalPeak(newLocalPeakidx);
localPeakValue = abs(localPeakValue(newLocalPeakidx));

	
