function [ timeData, yTimeData, frequencyData, amplitudeData, phaseData, approxCosineSeries, numberNonzeroTermsInApproxCosineSeries ] = FFTRealFunctionToCosines( amplitudeCutoffFactor )
%--------------------------------------------------------------
%   File:  FFTRealFunctionToCosines.m
% Authors: Carlos Osario and Paul Mitiguy
% Purpose: Periodic signal reconstruction using FFTCosine.
% License: Public domain.
%--------------------------------------------------------------

% Original signal parameters
% Note: Nyquist criteria says to sample at _HIGHER_ than 2x highest frequency.
highestExpectedFrequencyInHz   = 25;
desiredFrequencyResolutionInHz = 0.1;                              % Should be 5x-10x smaller than needed.
NyquistCriteriaSampleRateInHz = 4 * highestExpectedFrequencyInHz;  % 4x-16x is better than 2x.
timeStep = 1 / NyquistCriteriaSampleRateInHz; 
tFinal   = 1 / desiredFrequencyResolutionInHz;
numberOfElementsInSignal = NyquistCriteriaSampleRateInHz * tFinal;  
timeData = timeStep * ( 0 : numberOfElementsInSignal );   
frequencyData = NyquistCriteriaSampleRateInHz * linspace( 0, 1, numberOfElementsInSignal+1 );

% Create signal y(t).  To add noise, use:   0.6 * randn( size(timeData) );
% yTimeData = 3.0 * cos( 0.0 * 2*pi*timeData )  +  2.0 * cos( 1.0 * 2*pi*timeData ); 
% yTimeData = 2.0 * cos( 1.0 * 2*pi*timeData )  +  0.5 * cos( 8.0 * 2*pi*timeData ); 
% yTimeData = 0.5 * cos( 1.0 * 2*pi*timeData )  +  2.0 * cos( 8.0 * 2*pi*timeData ); 
yTimeData = 2.0 * cos( 2.0 * 2*pi*timeData )  +  1.0 * cos( 3.0 * 2*pi*timeData ); 
% yTimeData = 2.0 * cos( 2.75 * 2*pi*timeData )  +  1.0 * cos( 3.0 * 2*pi*timeData );   % Need desiredFrequencyResolutionInHz = 0.05;
% yTimeData = 1.0 * cos( 1.0 * 2*pi*timeData )  +  2.0 * cos( 2.0 * 2*pi*timeData )  +  1.0 * cos( 4.0 * 2*pi*timeData )  +  2.0 * cos( 8.0 * 2*pi*timeData ); 

% Fourier transform calculation - full spectrum.
Y = fft( yTimeData ) / numberOfElementsInSignal;

% fft returns its information in a "math" form, with each data point being a complex number.
% The phase of the complex number is the phase at that frequency.
% The amplitude of the complex number is the 1/2 amplitude at that frequency -- 
% unless the frequency is 0, in which case it is the full amplitude.
% To cull out the "engineering" useful amplitude and phase information, just
% view the first half of the data and scale all the amplitudes (other than the 
% 0-frequency data points) by a factor of 2.
halfWidthOfData = ceil( length(Y) / 2 );
Yhalf           = Y( 1 : halfWidthOfData );
phaseData       = angle( Yhalf );            % Phase of complex number          (complex number represents amplitude & phase).    
amplitudeData   = 2 * abs( Yhalf );          % Absolute value of complex number (complex number represents amplitude & phase).
if( frequencyData(1) == 0 ),  amplitudeData(1) = 0.5 * amplitudeData(1);  end;

% Choose frequencies of sufficient amplitude to reconstruct signal.
% Note:  0= < amplitudeCutoffFactor <= 1,  and is approximately 0.05
if( nargin==0 || length(amplitudeCutoffFactor)==0 || amplitudeCutoffFactor < 0.0 || amplitudeCutoffFactor > 1.0 ),  amplitudeCutoffFactor = 0.05;  end;
maxAmplitude = max( amplitudeData );
cutoffAmp = maxAmplitude * amplitudeCutoffFactor;
isYAmplitudeGreaterThanCutoffAmp = abs(Y) > cutoffAmp;  % Column matrix with 1 where condition is satisfied and 0 elsewhere.
Z = Y .* double( isYAmplitudeGreaterThanCutoffAmp );    % mask Y with zeros and ones

% Optional: Set cut-off frequency and implement low pass, high-pass, or band-pass filter.
% FcutoffLow  =   0;
% FcutoffHigh = 200;
% isFrequencyHigherThan = frequencyData > FcutoffLow  | HandleExtraJunkOnRhsOfY;
% isFrequencyLowerThan  = frequencyData < FcutoffHigh | HandleExtraJunkOnRhsOfY;
% isFrequencyInRange    = isFrequencyHigherThan  &  isFrequencyLowerThan; 
% Z = Z .* double(isFrequencyInRange);

% Return to the time domain - IFFT method
z = numberOfElementsInSignal * ifft( Z );

% Extract relevant frequencies and associated amplitude/phase to approximate signal with sum of cosine functions.
% Select single sided spectrum.
indexForTruncatedCosineSeries = amplitudeData > cutoffAmp;   % Column matrix with 1 where condition is satisfied and 0 elsewhere.
freqTruncatedCosines  = frequencyData( indexForTruncatedCosineSeries ); 
ampTruncatedCosines   = amplitudeData( indexForTruncatedCosineSeries );     
phaseTruncatedCosines = phaseData(     indexForTruncatedCosineSeries );     

% Return to the time domain - sum of cosines method ( approxCosineSeries should be nearly identical to z ).
numberNonzeroTermsInApproxCosineSeries = 0;
approxCosineSeries = zeros( 1, length(timeData) );
for i = 1 : length(freqTruncatedCosines);
   if( ampTruncatedCosines(i) > 0 ) 
      amplitudeMultiplier = 1.0;
      if( freqTruncatedCosines(i) == 0 ) amplitudeMultiplier = 1.0;  end;
      approxCosineSeries = approxCosineSeries  +  amplitudeMultiplier * ampTruncatedCosines(i) * cos( 2*pi*freqTruncatedCosines(i)*timeData + phaseTruncatedCosines(i) );
      numberNonzeroTermsInApproxCosineSeries = numberNonzeroTermsInApproxCosineSeries + 1;
   end
end

% Resize the frequencyData column matrix to half its size before returning.
frequencyData = frequencyData( 1: halfWidthOfData );

% Save information in text files.
timeDataToSave      = [ transpose(      timeData ),  transpose(     yTimeData ) ];
frequencyDataToSave = [ transpose( frequencyData ),  transpose( amplitudeData ),  transpose( phaseData ) ];
save( 'FFTRealFunctionToCosinesTimeData.txt',      'timeDataToSave',      '-ascii' );
save( 'FFTRealFunctionToCosinesFrequencyData.txt', 'frequencyDataToSave', '-ascii' );




