
%   
% An introduction to the basic steps involved in spectral analysis

% This routine performs a spectral analysis of a series of sine waves

% General scheme adapted from Muller and MacDonald, 2001, Ice Ages and
% Astronomical Causes

%% CREATE DATA

clear all;

t=[0:600];% creates a time vector  
%             starting at 0, going to 600.  the ' at the
%          end of the above statement converts a row vector to a column vector

%create 3 sine waves with various periods
r=(.2*sin(2*pi*t/10))+(.2*sin(2*pi*t/27))+(.2*sin(2*pi*t/50));

%detrend -- this tends to minimize long-term trends in the data (not so
%important in this case, but it is with real data sets)
r=detrend(r);
% next, plot the data just to look at it
figure(1);
plot(t,r); % plots the data in the time domain
%% FFT
  
% now we calculate the FFT after padding
npt=2^14; % this is the number of points, used below to pad the data vector with zeros
% so that the data set is a power of 2 in length -- makes the fast fourier transform routine runs faster
% this number should be about 10 times as large as the length of your data set

H=fft(r,npt);% this calls the fft routine that is built into Matlab and puts the result
% in a vector that contains the complex fourier coefficients, which consist
% of a real and imaginary part

% next we calculate what is called the Fourier power
P=H.*conj(H); % this gets rid of the imaginary part by multiplying every component of H by 
% its complex conjugate (real(H) - i*imag(H)) 

amp=abs(H); % sqrt(real(H).^2 + imag(H).^2) this is sometimes referred to as the Fourier amplitude
Pow=(amp.^2)./npt; % this is just another way of expressing the Power

%calculate frequency array -- assign frequencies to the Powers
Nyquist=abs(0.5/(t(2)-t(1))); % the Nyquist frequency is the highest frequency you can resolve given the data
f=linspace(0,Nyquist,npt/2);  % this creates a vector of the frequencies corresponding to the Power
% 

% normalize the Power -- this is just a convention to preserve the peak
% shape
P=P/mean(P(1:npt/2));
%plot
fmax=.2; % this reduces the frequencies that will be plotted
num=round(npt/2*fmax/Nyquist); % used to truncate the Power and freq


figure(2); % this creates a new figure window
plot(f(1:num),P(1:num)) % plots frequency on the x and Power on the y
xlabel('frequency')
ylabel('spectral power')
mark1x % this calles a routine that will print the period (1/f) 
% anywhere you click the mouse in the graph.

figure(3);
plot(f(1:num),amp(1:num)) % plots frequency on the x and Power on the y
xlabel('frequency')
ylabel('Amplitude')
mark1x

figure(4);
plot(f(1:num),Pow(1:num)) % plots frequency on the x and Power on the y
xlabel('frequency')
ylabel('Power Squared')
mark1x


%% FFT shortcut

[P,f]= dbfft(t,r,.2,2^14);

figure(3);
plot(f,P);
xlabel('frequency')
ylabel('spectral power')
mark1x 

%% Band-Pass Filtering

% isolating one compnent of a complex signal and then transforming it back
% into the time domain by means of an inverse FFT

%first define the frequency band (in terms of period)
hi=55; % i.e, the long period edge of the band
lo=45; % i.e, the short period edge of the band
   
f1=1/hi;
f2=1/lo;
n=length(t);
fc=fft(detrend(r));
fNyq=.5/abs(t(2)-t(1));

fre=linspace(0,2*fNyq,n)';

[m,k1]=min(abs(fre-f1)); % find k1 corresponding to f1
[m,k2]=min(abs(fre-f2)); % find k2 corresponding to f2
k3=n-k2+2; % conjugate points
k4=n-k1+2;
% zero all but the desired band
fc(1:k1)=zeros(k1,1);
fc(k2:k3)=zeros(k3-k2+1,1);
fc(k4:n)=zeros(k1-1,1);
%take inverse FFT; remove the imaginary part
rnew=real(ifft(fc));

figure(6);
plot(t,r,t,rnew);
legend('data','BP of data');

%% shortcut BP

rnew1=sqbandpass(t,r,55,45);

figure(7); plot(t,r,t,rnew1)
legend('data','BP of data');

%% Sample interval

t2=linspace (min(t),max(t),.1*length(t));
r2=interp1(t,r,t2);

figure(8);
plot(t,r,t2,r2);

[P2,f2]=dbfft(t2,r2,.04,2^14);

figure(9);
plot(f,P,f2,P2);
legend('data','down-sampled')
mark1x



