%function dOut=FFTDetrend(rawData)
%This function takes out a frequncy (e.g., seasonal period subtracted from 
%long-term monthly mean tide records) and the  harmonics (specify below as nHarm) 
%from a data file, by transforming it to frequency domain, subtracting the
%stated frequencies (as fSub below) and the nHarm harmonics on either side,
%then inverse transforming it.  The "de-seasoned" data are plotted atop the
%raw data. The data to process are in the second column of rawData.

dataIn=rawData(:,2);  %only the single column for fft
m=mean(dataIn);
cenData=dataIn-m; %centered data
fa=fft(cenData); %fft of raw data

fundamental=1/(length(cenData)/12); %in years if monthly data  %%%%%%MAYBE CHANGE%%%%%%%%%
nyq=1/(2/12); %also in years

freq=[fundamental:fundamental:nyq*2];

fSub=1;  %frequency to subtract.
nHarm=2; %number of Harmonics to subtract;
inyq=find(freq,nyq);
iFreqSub=[];
halfPeak=3; %half the peak width, so subtract this many points on eather side of target freq.
for i=1:nHarm+1 %subtract nHarm+the original freq
    temp = find(and(freq<i+halfPeak*fundamental,freq>i-halfPeak*fundamental));
    iFreqSub=[iFreqSub,temp,length(fa)-temp+2];
end

fb=fa;
fb(iFreqSub)=0;

subPlot(2,1,1)
plot(freq,abs(fa));
title('Spectrum overlaid with peaks removed')
hold all
plot(freq,abs(fb));
hold off

subplot(2,1,2);
dOut=ifft(fb)+m;
plot(cenData+m)
title('Raw Data, overlaid with Seasonally Detrended')
hold all
plot(dOut)
hold off