% EXAMPLE OF NOISY SINES WITH MEAN NORMALIZATION
a = 1.4;
b = 2;
cu = rand(1,N);
y2 = a*sin(b*t)+cu*amp +1; %Unif error
a2 = 0.4;
b2 = 2;
amp = 2.3;
y3 = a2*sin(b2*t)+cu(randperm(N))*amp - 2; %Unif error
plot([1:N], y2);
hold on
plot([1:N], y3);

% EXAMPLE OF AUTOCORRELATION COMPUTATION

% Generate noisy sine
N = 100;
t = (0:4*pi/N:4*pi);
N = N+1;
a = 3;
b = 2;
cn = randn(1,N);
amp = 2;

% generate data and calculate "y"
y1 = a*sin(b*t)+cn*amp; %Gaussian error

% plot results
plot(t, y1);

% plot cross covariance
[cov_ww,lags] = xcov(y1);
stem(lags,cov_ww)


%%%%%%%% FOURIER ANALYSIS

Fs = 150; % Sampling frequency
t = 0:1/Fs:1; % Time vector of 1 second 
f = 5; %sine wave of f Hz

%SQUARE SIGNAL
%x = square(2*pi*t*f); 

%Noisy sines of 15 and 40 Hz
%x = sin(2*pi*15*t) + sin(2*pi*40*t); 

%Cosine wave
%pha = 1/3*pi; %phase shift
%x = cos(2*pi*t*f + pha); 

%Chirp signal
x = chirp(t,0,1,Fs/6); 

nfft = length(x); % Length of FFT
% Take fft, padding with zeros so that length(X) is equal to nfft 
X = fft(x,nfft); 
%FFT is symmetric, you could throw away second half 
%X = X(1:nfft/2); 
% Take the magnitude of fft of x
mx = abs(X);
% Frequency vector
f = (0:nfft-1)*Fs/nfft; %/2
% Generate the plot, title and labels. 
figure(1);
plot(t,x);
title('Signal'); 
xlabel('Time (s)'); 
ylabel('Amplitude');
figure(2);
plot(f,mx); 
title('Power Spectrum'); 
xlabel('Frequency (Hz)'); 
ylabel('Power'); 

%Frequency filtering
Xlp = X; 
Xlp(f>=30) = 0; %Filter out higher frequencies than 30Hz
ylp = ifft(Xlp,'symmetric');
plot(t,ylp);
title('Signal+Fourier Reconstruction'); 
xlabel('Time (s)'); 
ylabel('Amplitude');