
%% 24/01/2018 -- christophe@pallier.org

%% plot a sine wave (sinusoide)

freq = 0.3
phase= 0
amplitude = 3

time_to_plot = 10 % seconds
sample_rate = 100 % samples per second
num_samples = sample_rate * time_to_plot

time = 0:(1/sample_rate):time_to_plot;

signal = amplitude*sin((2*pi*freq*time) - phase);

plot(time, signal)

%% plot two sinusoides, varying amplitude, freq, and phase

freq1 = 1; amp1 = 2; phase1 = 0;
freq2 = 1; amp2 = 0.5; phase2 = 0.5;

y1 = amp1*sin((2*pi*freq1*time) - phase1);
y2 = amp2*sin((2*pi*freq2*time) - phase2);
figure()
plot(time, y1, 'r', time, y2, 'b')

%%

amp1=3; freq1=2; phase1=0;
amp2=1; freq2=0.5; phase2=0;
y1 = amp1*sin((2*pi*freq1*time) - phase1);
y2 = amp2*sin((2*pi*freq2*time) - phase2);
figure()
plot(time, y1, 'r', time, y2, 'b')

%% 
amp1=1; freq1=2; phase1=0;
amp2=1; freq2=2; phase2=pi/2;
y1 = amp1*sin((2*pi*freq1*time) - phase1);
y2 = amp2*sin((2*pi*freq2*time) - phase2);
figure()
plot(time, y1, 'r', time, y2, 'b')

%% summing sine waves

amp1=3; freq1=12; phase1=0;
amp2=3; freq2=11; phase2=0;
y1 = amp1*sin((2*pi*freq1*time) - phase1);
y2 = amp2*sin((2*pi*freq2*time) - phase2);
figure()
plot(time, y1 + y2)

%%

amp1=3; freq1=2; phase1=0;
amp2=1; freq2=1; phase2=pi/4;
y1 = amp1*sin((2*pi*freq1*time) - phase1);
y2 = amp2*sin((2*pi*freq2*time) - phase2);
figure()
plot(time, y1 + y2)


%% summing more than 2 sine waves (6 here)

f = [1,2,4,8,16,32]
a = [1,0.5,0.8,2,0.3,0.5]
p = [0,1,0,1,0,0]

y = zeros(1, size(time, 2));
for i = 1:size(f, 2)
    y = y + a(i) * sin(2*f(i)*time-p(i));
endfor

plot(time, y)

%% Decomposition


% basis functions

time = 0:0.005:5;
s0 = ones(1, size(time, 2));
s1 = sin(2*pi*time); c1 = cos(2*pi*time);
s2 = sin(2*pi*2*time); c2 = cos(2*pi*2*time);
s3 = sin(2*pi*4*time); c3 = cos(2*pi*4*time);
s4 = sin(2*pi*8*time); c4 = cos(2*pi*8*time);
s5 = sin(2*pi*16*time); c5 = cos(2*pi*16*time);
s6 = sin(2*pi*32*time); c6 = cos(2*pi*32*time);
s7 = sin(2*pi*64*time); c7 = cos(2*pi*64*time);

X = [s0; s1; c1; s2; c2; s3; c3; s4; c4; s5; c5; s6; c6; s7; c7]';

plot(time, X(:,1:5))

image(X*100)

%%  Decomposing a signal on the basis functions using multiple linear regression

% let us create a saw tawtooth function
pkg load signal  % hum, works for octave, but likely not in matlab!
%% see http://www.matrixlab-examples.com/periodic-function.html
y = sawtooth(2 * pi * time);
plot(time, y)

betas = X \ y';

reconstructed_signal = X * betas;
plot(time, reconstructed_signal)

%% foutier tranform (perdiodogram)
 
plot(Pxx)

%% Spectrogram 

dt = 0.0005
t = 0:.0005:20;
s1 = sin(2*pi*100*t);
s2 = 2*sin(2*pi*400*t);

% create a transient "chirp"
mask = t > 10 & t < 12;
s2 = s2 .* mask;

% add some noise into the mix
%nse = 0.01*rand(size(t,2))

x = s1 + s2;   % the signal
NFFT = 1024       % the length of the windowing segments
Fs = floor(1.0/dt)  % the sampling frequency

plot(t, x)
figure()
specgram(x, NFFT=NFFT, Fs=Fs)
