% (1) Generere et datasett med n målinger, hvor målingene varierer sinusoidalt med amplitude A. Det skal være en normalfordelt støy med standardavvik s=0.2*A, og en liten linær drift i målingene fra første til siste måling.
% Set up n=1000 samples taken at sampling frequency f_s =1\text{ }kHz of a sinus signal with amplitude A=2\ldotp 5 V and choose a frequency that will give us N_p =6 periods within the sampling time of \frac{n}{f_s }=1 s.

fs=1000; % Hz
n=1000;
A=2.5;
N_p=6;
t_total = n/fs
f=6/t_total

%
% To this signal we will add gaussian random noise with an amplitude that is 20% the size of the signal amplitude and a linear drift over the whole measurement time that is twice the amplitude.

noise_sigma = 0.2*A
kdrift=2*A
noise=noise_sigma*normrnd(0,1,1,n); % normrnd(mu,sigma,M,N)
% alternatively stoy = normrnd(0,sigma_stoy*A,1,n);

%
% (2) Plot signalet. Benevning på y-aksen skal være Volt [V]. Langs x-aksen skal benevningen være sekunder, og signalet er tatt opp med en frekvens på 1kHz.


t=linspace(0,t_total,n);
drift=kdrift*t/t_total;

sig = drift + noise + A*sin(2*pi*f*t);

plot(t,sig,'.')
title('(2) Measured voltage versus time')
xlabel('tid [s]')
ylabel('Volt [V]')

%
% (3) Lag et histogram av målingene i en ny figur med passe antall kolonner i forhold til antallet målinger n. Husk benevning på aksene. Kommenter på fordelingen.
% I played around with the number of samples and bins and this simple parameterization gave a reasonable comprimise between seeing the details of the shape of the distribution versus not having too large bin-to-bin fluctuations.

nbins=round(8*log10(n))
histogram(sig,nbins);
title('(3) Histogram of sampled voltage')
xlabel('Sampled voltage [v]')
ylabel('Number of samples')
% The distribution of sampled voltages doesn't look much like a normal distribution.

%
% (4) Lag et histogram som viser fordelingen av kun støyen i signalet, dvs signalet trukket fra sinus- og drift-delen.
% Since we have stored the noise in a separate variable we can simply histogram it.

histogram(noise,nbins);
title('(4) Histogram of noise')
xlabel('Pure noise component [V]')
ylabel('Number of samples')

%
% (5) Normaliser arealet til kolonnene (dvs at det samlede arealet til kolonnene skal være 1), og plot i en ny figur (figure 4) punktene  som representerer den normaliserte høyden på kolonnene sammen med normalfordelingskurven med standardardavvik s. Lek litt med forskjellige verdier av n for å se hvor mange målinger som skal til før vi har en god overenstemmelse med normalfordelingskurven.
% We can do the normalization by hand or we can use the properties of the stored noise histogram. I chose to use the properties of the histogram. I extract the centers of the bins and the height of the pdf-normalized histogram and overlay the normal pdf with the known parameters of the noise.

% Create a new version of the noise histogram and change its normalization to pdf
hn=histogram(noise,nbins)
hn.Normalization='pdf';
title('Ignore this temporary histogram')

% Extract the centers and heights of the bins a
% 
x = hn.BinEdges(1:end-1) + hn.BinWidth/2;
y = hn.Values
binwidth = hn.BinWidth;

% Calculate the normal pdf on a finer grid
xn = linspace(hn.BinEdges(1),hn.BinEdges(end),100);
yn = normpdf(xn,0,noise_sigma);

% Plot the sampled noise and overlay the normal pdf
plot(x,y,'k.','MarkerSize',25)
hold on
plot(xn,yn,'r-')
title('(5) Normalisert stoy og norm-pdf')
xlabel('Sampled voltage [V]')
ylabel('Probability density')
legend('Sampled noise','Predicted noise')
hold off

%
% (6) Se om du kan finne en funksjon som finner den linære driften i målingene, og se om du klarer å bruke dette til å manipulere bort denne fra målingene. Her holder det å gi en kommentar med forslag til fremgangsmåte.
% We shall use the Matlab function detrend, which fits the data to a straight line and then subtracts this line from the data. Note that the drift is relatively small and the underlying signal is periodic and there are few periods and the number of periods is not an integer the subtraction will not work as well as it would for an integer number of periods or if the number of periods is large.

sig_detrended = detrend(sig,'linear');
plot(t,sig_detrended,'k.')
title('(6) Measured voltage versus time, drift-corrected')
xlabel('tid [s]')
ylabel('Volt [V]')

% Overlay the predicted signal
hold on
plot(t,A*sin(2*pi*f*t),'r-')
legend('Corrected samples','Predicted signal')