######################################################
## Intro to R and Fourier Analysis 
##
##1. Basic data structure
##2. For loop
##3. Functions for basic operations
##4. Plots and graphs
##
##
#######################################################

## Create a cosine vector of length=Tlen and frequency=freq
Tlen = 128; freq=2;
timet = seq(0:(Tlen-1));
cos2 = cos(2*pi*freq*timet/Tlen);

## Plot this vector
plot( x=timet, y=cos2, type="l");

## Create a set of sine and cosine vectors
Trigmatrix = matrix(ncol=Tlen, nrow=Tlen);

## The cosine waves occupy columns 1 to (Tlen/2 + 1)
for (freq in (0:(Tlen/2))){
colindex = freq + 1;
Trigmatrix[, colindex] = cos(2*pi*freq*timet/Tlen);
} 

## The sine waves occupy columns (Tlen/2 + 1) to T.
for (freq in (1:(Tlen/2 -1))){
colindex = Tlen/2 + 1 + freq;
Trigmatrix[, colindex] = sin(2*pi*freq*timet/Tlen);
} 

## Plot some of the cosine and sine waves
par(mfrow=c(2,3)); 
ts.plot(Trigmatrix[, 1]); title(paste("Cos 0")); 
ts.plot(Trigmatrix[, 2]); title(paste("Cos 1")); 
ts.plot(Trigmatrix[, 10]); title(paste("Cos 9"));
ts.plot(Trigmatrix[, Tlen/2 + 30]); title(paste("Sine 29")); 
ts.plot(Trigmatrix[, (Tlen/2 + 2)]); title(paste("Sine 1")); 
ts.plot(Trigmatrix[, (Tlen/2 + 10)]); title(paste("Sine 9")); 

## Are the cosine and sine waves orthogonal? 
## Check if (Trigmatrix)' (Trigmatrix) is a diagonal matrix

IP = t(Trigmatrix)%*%Trigmatrix;

## Then examine the off diagonals of IP. One way is to plot 
## each row. Check if only the r-th element of row r is 
## large and the rest are essentially zero.

par(mfrow=c(2,2)); 
ts.plot(IP[, 1]); title(paste("< Cos 0, ALL >"));
ts.plot(IP[, 2]); title(paste("< Cos 1, ALL >"));
ts.plot(IP[, (Tlen/2 + 2)]); title(paste("< Sin 1, ALL >"));
ts.plot(IP[, (Tlen)]); title(paste("< Sin 63, ALL >"));

## Other commands
## Inner product between Cos0 and Cos1 waveforms
Cos0 = Trigmatrix[,1];
Cos1 = Trigmatrix[,2];
sum(Cos0*Cos1)

## Generate a time series using arima.sim
## AR(1) with parameter 0.9
x = arima.sim(Tlen, model=list(ar=c(0.9), ma=c(0)));
par(mfrow=c(1,1));
ts.plot(x);

## For plotting purposes, we'll rescale x
xnew = x - mean(x);
xnew = 2*xnew/((max(x) - min(x)));
ts.plot(xnew)

## Plot xnew with the sine and cosine waves

par(mfrow=c(2,2));
Cos3 = Trigmatrix[,4]; plot(xnew, ylim=c(-1,1)); lines(Cos3, lty=3); 
title(paste("Cos Low Freq"));
Sin3 = Trigmatrix[,(Tlen/2 + 4)]; plot(xnew); lines(Sin3, lty=3); 
title(paste("Sin Low Freq")); 

Cos45 = Trigmatrix[,46]; plot(xnew); lines(Cos45, lty=3); 
title(paste("Cos High Freq"));
Sin45 = Trigmatrix[,(Tlen/2 + 47)]; plot(xnew); lines(Sin45, lty=3); 
title(paste("Sin High Freq")); 


x = arima.sim(Tlen, model=list(ar=c(-0.9), ma=c(0)));
## Compute the coefficients corresponding to the Trig waveforms 
Coscoef = c(1:(Tlen/2+1));
for (k in 1:(Tlen/2 + 1)){
Coswave = Trigmatrix[,k];
Normwave = sum(Coswave*Coswave);
Coscoef[k] = (1/Normwave)*sum(Coswave*x);
}

Sincoef = c(1:(Tlen/2 -1));
for (k in 1:(Tlen/2-1)){
Sinwave = Trigmatrix[,(Tlen/2 + 1 +k)];
Normwave = sum(Sinwave*Sinwave);
Sincoef[k] = (1/Normwave)*sum(Coswave*x);
}

## Sum of squares of Coefficients
Amp = c(1:(Tlen/2+1));
Amp[1] = (Coscoef[1])^2;
Amp[(Tlen/2+1)] = (Coscoef[(Tlen/2+1)])^2;
Amp[2:(Tlen/2)] = (Coscoef[2:(Tlen/2)])^2 + (Sincoef[1:(Tlen/2 -1)])^2;

## Plot the Sum of Squares of coefficients

par(mfrow=c(1,1)); 
plot(Amp, type="p");