function FTdemoFeb14a

% For demonstrasjon i FYS2130 7. februar 2014, Arnt Inge Vistnes
% Ved bruk: Kommenter ut en "return" etter en annen etter som vi skrider
% fram i bruk av programmet. Kommenter ut hele blokker når det er aktuelt.
clear all;

% Leser inn et lydsignal, plotter ene kanalen og spiller lyden
N = 2^11;
fs = 44100;
T = N/fs;
t = linspace(0,T,N);
s = 'piclDet.wav';
[ffull, Fs] = audioread(s,[1 N]);
f = ffull(:,1);  % Henter ut venstre (eller høyre) kanal only
plot(t,f,'-b');
xlabel('Tid (s)')
ylabel('Lydsignal (rel verdi)');
sound(ffull,Fs);

return;

% Leser av frekvens vi vil analysere med, og plotter f,g og f*g slik som
% tidligere i de første figurene i kapittel 4 (Kommenter ut etter bruk!)
plot(t,f,'-b');
title('Funksjon blå, exp(i omega t) rød');
freq = 602.9;  % 28/0.04644 lest av tidsbildet av signalet
g = cos(2*pi*freq*t);   % Velg cos for realdel og -sin for imagdel
hold on;
plot(t,g,'-r');

figure,
plot(t,f.*g','-k');
title('Produktet av f og exp(i omega t)');

Ca_integrert = sum(f.*g')   % Skrives til skjerm

return;

% Foretar en FFT og plotter realdel/imaginærdel pluss abs-verdi
% Merk symmetriene, realverdien symmetrisk rundt halve samplingsfrekvensen,
% mens imaginærdelen er antisymmetrisk
F = fft(f);
figure;
freq = linspace(0,fs*(N-1)/N, N);
plot(freq,real(F),'-r');
hold on;
plot(freq,imag(F),'-b');
xlabel('Frekvens (Hz)')
ylabel('Fourierkoeffisienter (rel verdier)');
legend('real(F): ', 'imag(F): ');
figure;
plot(freq,abs(F),'-k');xlabel('Frekvens (Hz)')
ylabel('Fourierkoeffisient (rel verdier)');
title('Absloluttverdien av frekvensspekteret');

return;

% Invers fouriertransform av opprinnelig fourierspekter
H = ifft(F);
figure;
Hv = H(:,1);
plot(Hv);

return;

% Leser av koeffisienter fra FT og forsøker å etterligne opprinnelig signal
% litt etter litt (alt etter hvor mange fourierkoeffisienter vi velger å ta
% med)
figure;
freq = 602.9;
signal = -334.5*cos(2*pi*freq*t) - 333.5*sin(2*pi*freq*t);
% freq = 1206;
% signal = signal + 61.67*cos(2*pi*freq*t) - 36.33*sin(2*pi*freq*t);
% freq = 1809;
% signal = signal - 29.9*cos(2*pi*freq*t) + 124.6*sin(2*pi*freq*t);
% freq = 2412;
% signal = signal - 59.78*cos(2*pi*freq*t) + 15.17*sin(2*pi*freq*t);
% freq = 2993;
% signal = signal - 0*cos(2*pi*freq*t) + 11.65*sin(2*pi*freq*t);
% freq = 64.6;
% signal = signal - 0*cos(2*pi*freq*t) -  24.61*sin(2*pi*freq*t);
% freq = 21.53;
% signal = signal - 0*cos(2*pi*freq*t) + 27.09*sin(2*pi*freq*t);
signal = 2*signal/N;   % Normering må sjekkes manuelt. Ulikt i ulike prog.språk
plot(t,f,'-b');
xlabel('Tid (s)')
ylabel('Lydsignal (rel verdi)');
hold on;
plot(t,signal,'-r');
legend('Oppr. signal', 'Etterlignet signal');

return;

% En liten morsomhet til slutt: Hvorfor ikke spille den lydsnutten vi nå har
% generert?
sound(signal,fs);

