function varargout = RECSINEWAVES( Data, Time, number )
%   RECSINEWAVES Recover the sine waves of the waveform
%   
%      [X,Y] = RECSINEWAVES(DATA,TIME, NUMBEROFTONES);
%   
%    This function will find the local peaks in the fft of the data assuming the 
%    highest peak is the first peak. Then rebuild a waveform using only those frequencies.
%   
%      DATA is the data from the instrument.
%  
%      TIME is a vector of time points where data was taken.
%   
%      NUMBEROFTONES is the maximum number of tones to use in reconstructing
%      the waveform.  If the maximum number is found is less than this then
%      only the tones found are used.  If NUMBEROFTONES is not provided then the 
%      peaks in the spectrum above the noise floor are used.
%   
%    Example:
%      [x,y] = RECSINEWAVES(Data,Time);
%   
%    See also
%    RECPRIMARYWAVE, REMPRIMARYWAVE 


if nargout ==3
    varargout{1} = 'Time [sec]';
    varargout{2} = 'Amplitude [V]';
    varargout{3} = 'Recover the sine waves of the waveform';
    return;
end;

if nargin==0
    help(mfilename)
    return;
end;


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

if isempty(val)
    fData = zeros(size(Data));
else
    numTones = length(val);
    if nargin==3
        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);
    freq = freq(:)';
    fData=zeros(size(freq));
    
    idx=min(length(val),20); % This is just to make it faster (DEMO Code).
    
    for c=1:idx;
        fData=fData+val(c)*sin(freq*freqidx(c));
    end;
end;
varargout{1} = Time;
varargout{2} = fData;

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));



Return to Content