%*********************************************************************
%*********************************************************************
%********* M1 " DSV mit Matlab - Anwendung der DFT/FFT" **************
%*********************************************************************
%*********************************************************************

echo on;
%*********************************************************************
% 1 Approximation der Fourierreihe eines periodischen Rechtecksignales
%   mittels FFT
%*********************************************************************


% --- 1.1 Vorbereitung
%*********************************************************************
echo off;
% Berechnen Sie zunächst den Betrag der komplexen Fourierkoeffizienten des
% in Bild 1 skizzierten periodischen Rechtecksignals mit normierter
% Amplitude A=1000 und tragen Sie den Wert der ersten 
% 11 Fourierkoeffizienten |ck| für k=0,1...10 in Tabelle 1 ein.

A = 1000;
k=[0:15];
% deltaT/Tp = 1/4
% Yk = A * deltaT/Tp * si(pi*k*deltaT/Tp)
Yk = A * 1/4 * sinc(k*1/4);
Yk = abs(Yk)'; k = k';

% Tabelle Ausgeben
f = figure;
data = [k Yk];
colnames = {'k', '|Yk|'};
t = uitable(f, 'Data', data, 'ColumnName', colnames, ...
                   'Position', [20 100 185 310]);
%pause;
close all;

echo on;

%% --- 1.2 des diskreten Mustersignals y1 per Simulation
%*********************************************************************
echo off;
% Es soll ein Zeitausschnitt über eine Periode der Dauer TP (siehe Bild 1)
% diskret mit L=16 Samples erzeugt werden:
% ==> Welche Abtastfrequenz ist hierfür notwendig:
L = 16;
Tp = 10^-3; % -> TP = 0 bis 10^-3 -> f = 1000 Hz
echo on;
% fa1 = L/Tp = 16/10^-3 = 16 kHz
echo off;
fa1 = 16000;
L=16;
k=[0:L-1]; %Index
t=k*1/fa1; % Zeitvektor d. Länge 1 ms
% Ausgangsvektor y, erzeugt durch logische
% Verknüpfung mit dem Zeitvektor t:
y=[t<0.25e-3];
% heißt: es wird ein y Vektor mit gleicher
% Länge wie t erzeugt und für jedes Element
% in t die Bedingung t<0.25e-3 geprüft; an den
% Indizes, wo Bedingung erfüllt ist, wird y=1; sonst 0
y=A*y; %Amplitude berücksichtigen
stem(t,y); xlabel('t/ms \rightarrow'); ylabel('y_k');
title('1.2 Zeitdiagramm des abgetasteten Rechteckimpulses',... 
  'FontWeight','bold')
pause;
close all;

echo on;

%% --- 1.3 Analyse des Mustersignals y1 (eine Periode) mit (DFT) FFT
%*********************************************************************
echo off;

% Da eine Länge von L=16=2^4 für das Mustersignal gewählt wurde, 
% wird hier der FFT Algorithmus verwendet.
% Es sollen die in Tabelle 1 schon theoretisch ermittelten 
% Fourierkoeffizienten nun durch FFT approximiert werden:

korr = 1/L 
Yd1=korr*abs(fft(y)) % korr muss noch bestimmt werden
f=[0:L-1]*fa1/L

stem(f,Yd1) % Plot mit diskreten Linien
xlabel('f/Hz \rightarrow'); ylabel('Yd1');
axis([0 fa1 0 250]);grid;
title('1.3 Linienspektrum des abgetasteten Rechteckimpulses',... 
  'FontWeight','bold')

% Tabelle Ausgeben
figure(2);
f = figure;
Yd1 = abs(Yd1)';
k = k';
data = [k Yd1];
colnames = {'k', '|Yd1|'};
t = uitable(f, 'Data', data, 'ColumnName', colnames, ...
                   'Position', [20 100 185 310]);
pause;
close all;

echo on;
%% --- 1.4 Veränderung der Abtastfrequenz: Mustersignal y2 (1 Periode)
%*********************************************************************

% Mit welcher Maßnahme lassen sich eben festgestellte Abweichungen
% von |Yd1| reduzieren? -> Abtastfrequenz erhöhen!
echo off;
fa2 = fa1*4
L2 = L*4


k=[0:L2-1];
t=k*1/fa2;
y=[t<0.25e-3];
y=A*y; %Amplitude berücksichtigen

% Zeitdiagramm
stem(t,y); xlabel('t/ms \rightarrow'); ylabel('y_k');
title('1.4 Zeitdiagramm des abgetasteten Rechteckimpulses bei erhöhtem fa und L',... 
  'FontWeight','bold')

% Linienspektrum
Yd2=1/L*abs(fft(y));
f=[0:L2-1]*fa2/L2;
figure;
stem(f,Yd2); xlabel('f/Hz \rightarrow'); ylabel('Yd2');
title('1.4 Linienspektrum des abgetasteten Rechteckimpulses bei erhöhtem fa und L',... 
  'FontWeight','bold')

% Tabelle Ausgeben
figure;
f = figure;
Yd2 = abs(Yd2)';
k = k';
data = [k Yd2];
colnames = {'k', '|Yd2|'};
t = uitable(f, 'Data', data, 'ColumnName', colnames, ...
                   'Position', [20 100 210 310]);
pause;
close all;

echo on;
%% --- 1.5 Verlängerte Messzeit über 4 Perioden : Mustersignal y3
%*********************************************************************
echo off;
% Was bringt es, wenn die ursprüngliche Folge aus Abb. 1 statt mit einer 
% Periode beispielsweise mit 4 Perioden dargestellt wird?
% (Abtastfrequenz fa wie ganz am Anfang!!) -> höhere grafische Auflösung

fa1= 16000;
L=16;
A=1000;
k=[0:L-1]; 
t=k*1/fa1;
y=[t<0.25e-3];
y=A*y;

% 4* hintereinander gesetzt; Länge jetzt 4*L
y3=[y y y y];
L=4*L;

% Zeitdiagramm
k=[0:L-1];
t=k*1/fa1;
stem(t,y3); xlabel('t/ms \rightarrow'); ylabel('y_k');
title('1.5 Zeitdiagramm des abgetasteten Rechteckimpulses (4 Perioden)',... 
  'FontWeight','bold')

% Linienspektrum
figure;
Yd1=1/L*abs(fft(y3));
f=[0:L-1]*fa1/L;
stem(f,Yd1); xlabel('f/Hz \rightarrow'); ylabel('Yd2');
title('1.5 Linienspektrum des abgetasteten Rechteckimpulses (4 Perioden)',... 
  'FontWeight','bold')

% Welche Auswirkungen ergeben sich im Frequenzbereich ?
% - allgemein: Mehr Abtastpunkte
% - für die approximierten Fourierkoeffizienten: nicht genauer, da fa nicht
%   erhöht wurde

pause;
close all;

echo on;
%% --- 1.6 Auswertung einer nicht-ganzzahligen Anzahl von Perioden/ y4
%*********************************************************************
echo off;

L=length(y3);
L=L-4;
y4=y3(1:L);

% Verändert sich etwas an der Energie des Signals? -> Nein

% Zeitspektrum
k=[0:L-1];
t=k*1/fa1;
stem(t,y4);

% Linienspektrum
figure;
Yd4=1/L*abs(fft(y4));
f=[0:L-1]*fa1/L;
stem(f,Yd4); xlabel('f/Hz \rightarrow'); ylabel('Yd2');
title('1.6 Linienspektrum des abgetasteten Rechteckimpulses (Abgeschnitten)',... 
  'FontWeight','bold')

% Welche Veränderung ergibt sich im Linienspektrum?
% -> Leckeffekt + picket fence

pause;
close all;

echo on;
%*********************************************************************
% 2 Fourier-Spektrum eines einzelnen Rechteckimpulses
%*********************************************************************


%% --- 2.1 Vorbereitung
%*********************************************************************

% Y(0) = 250 uV/Hz
% Y(0)* deltaT = 250 

% --- 2.2 Approximation des kontinuierlichen Amplitudenspektrums mittels FFT
%*********************************************************************
echo off;
% Was muss man tun, damit man nach der FFT in der graphischen Darstellung 
% des approximierten kontinuierlichen Spektrums eine gute Auflösung erzielt ?
% -> Messdauer N*Ta muss mindestens T0 sein

% Wie muss das FFT-Ergebnis skaliert werden, um eine physikalisch 
% korrekte Approximation zu erreichen: -> x[k]/fa

close all;
fa1= 64000; % ?ausfüllen!
L=2^9;
A=1000;
k=[0:L-1]; %Index
t=k*1/fa1; % Zeitvektor d. Länge 1 ms
% Ausgangsvektor y, erzeugt durch logische
% Verknüpfung mit dem Zeitvektor t:
y=[t<0.25e-3];
% heißt: es wird ein y Vektor mit gleicher
% Länge wie t erzeugt und für jedes Element
% in t die Bedingung t<0.25e-3 geprüft; an den
% Indizes, wo Bedingung erfüllt ist, wird y=1; sonst 0
y=A*y; %Amplitude berücksichtigen
stem(t,y); xlabel('t/us \rightarrow'); ylabel('y_k');

figure;
Yd1=1/fa1*abs(fft(y));
f=[0:L-1]*fa1/L;
plot(f(1:256),Yd1(1:256));
hold 
Yk=abs(A*250e-6*sinc(f*250e-6));
plot(f(1:256),Yk(1:256),'ro-'); xlabel('f/Hz \rightarrow');
legend('Approcimation', 'Exakt', 1);

pause;
close all;

echo on;
%*********************************************************************
% 3 Detektion von Sinussignalen, Fensterung
%*********************************************************************


%% --- 3.1 Vorbereitung
%*********************************************************************

% N = 256 -> fa = 2560

% Was für ein Signal wird hier erzeugt? Frequenzen? Amplituden?
% Sinusschwungung + Rauschen; Amplituden: 100 und 2; 

% Welcher Pegelunterschied der zwei Sinusanteile in dB ergibt sich?
% -> 20log(A1/A2) = 34dB

% --- 3.2 Fall A; Abtastfrequenz ist frei wählbar
%*********************************************************************
echo off;

N=256;
fa=2560;
k=[0:N-1];
y=100*sin(pi*1980*k/fa)+2*sin(pi*2020*k/fa)+2*rand(1,N)-1;
Yk=abs(fft(y));
Ymax=max(Yk) %Maximum feststellen
Yk=20*log10(Yk/Ymax);
f=k*fa/N;
plot(f,Yk,'bo-'); grid; axis([0 fa/2 -100 0]); xlabel('f/Hz \rightarrow');
ylabel('Yk');
title('3.2 Lienienspektrum der Sinusschwingung + Rauschen',... 
  'FontWeight','bold')

pause;
close all;

echo on;
%% --- 3.3 Fall B; Abtastfrequenz ist fest vorgegeben
%*********************************************************************
echo off;
fa = 4000;
N=400;
k=[0:N-1];
y=100*sin(pi*1980*k/fa)+2*sin(pi*2020*k/fa)+2*rand(1,N)-1;
Yk=abs(fft(y));
Ymax=max(Yk); %Maximum feststellen
Yk=20*log10(Yk/Ymax);
f=k*fa/N;
plot(f,Yk,'bo-'); grid; axis([0 fa/2 -100 0]);
xlabel('f/Hz \rightarrow'); ylabel('Yk');
title('3.3 Linienspektrum der Sinusschwingung + Rauschen (fa =4000Hz)',... 
  'FontWeight','bold')
hold;

N=100;
k=[0:N-1];
y=100*sin(pi*1980*k/fa)+2*sin(pi*2020*k/fa)+2*rand(1,N)-1;
Yk=abs(fft(y));
Ymax=max(Yk); %Maximum feststellen
Yk=20*log10(Yk/Ymax);
f=k*fa/N;
plot(f,Yk,'ro-');  legend('N = 400', 'N = 100', 1);

pause;
close all;

fa=4000;
N=256;
k=[0:N-1];
y=100*sin(pi*1980*k/fa)+2*sin(pi*2020*k/fa)+2*rand(1,N)-1;
y=y.*hann(N)';
Yk=abs(fft(y));
Ymax=max(Yk) %Maximum feststellen
Yk=20*log10(Yk/Ymax);
f=k*fa/N;
plot(f,Yk,'ro-'); grid; axis([0 fa/2 -100 0]);
xlabel('f/Hz \rightarrow'); ylabel('Yk');
title('3.3 Linienspektrum fa = 256 mit HANN-Fenster',... 
  'FontWeight','bold')
hold;

N=256;
k=[0:N-1];
y=100*sin(pi*1980*k/fa)+2*sin(pi*2020*k/fa)+2*rand(1,N)-1;
Yk=abs(fft(y));
Ymax=max(Yk); %Maximum feststellen
Yk=20*log10(Yk/Ymax);
f=k*fa/N;
plot(f,Yk,'bo-'); legend('HANN', 'Original', 1);

pause;
close all;

N=512;
k=[0:N-1];
y=100*sin(pi*1980*k/fa)+2*sin(pi*2020*k/fa)+2*rand(1,N)-1;
Yk=abs(fft(y));
Ymax=max(Yk); %Maximum feststellen
Yk=20*log10(Yk/Ymax);
f=k*fa/N;
plot(f,Yk,'bo-');
hold

y=y.*hann(N)';
Yk=abs(fft(y));
Ymax=max(Yk); %Maximum feststellen
Yk=20*log10(Yk/Ymax);
f=k*fa/N;
plot(f,Yk,'ro-'); grid; axis([0 fa/2 -100 0]);
xlabel('f/Hz \rightarrow'); ylabel('Yk');
title('3.3 Linienspektrum fa = 512 mit HANN-Fenster',... 
  'FontWeight','bold');  legend('HANN', 'Original', 1);

pause;
close all;
figure;
N=512;
k=[0:N-1];
y=100*sin(pi*1980*k/fa)+2*sin(pi*2020*k/fa)+2*rand(1,N)-1;
y=y.*bartlett(N)';
Yk=abs(fft(y));
Ymax=max(Yk); %Maximum feststellen
Yk=20*log10(Yk/Ymax);
f=k*fa/N;
plot(f,Yk,'bo-'); grid; axis([0 fa/2 -100 0]);
xlabel('f/Hz \rightarrow'); ylabel('Yk');
title('3.3 Linienspektrum mit BARTLETT-Fenster',... 
  'FontWeight','bold'); 

pause;
close all;
figure;
N=512;
k=[0:N-1];
y=100*sin(pi*1980*k/fa)+2*sin(pi*2020*k/fa)+2*rand(1,N)-1;
y=y.*bartlett(N)';
Yk=abs(fft(y));
Ymax=max(Yk); %Maximum feststellen
Yk=20*log10(Yk/Ymax);
f=k*fa/N;
plot(f,Yk,'bo-'); grid; axis([0 fa/2 -100 0]);
xlabel('f/Hz \rightarrow'); ylabel('Yk');
title('3.3 Linienspektrum mit BARTLETT-Fenster',... 
  'FontWeight','bold'); 

pause;
close all;
figure;
N=512;
k=[0:N-1];
y=100*sin(pi*1980*k/fa)+2*sin(pi*2020*k/fa)+2*rand(1,N)-1;
y=y.*blackman(N)';
Yk=abs(fft(y));
Ymax=max(Yk); %Maximum feststellen
Yk=20*log10(Yk/Ymax);
f=k*fa/N;
plot(f,Yk,'bo-'); grid; axis([0 fa/2 -100 0]);
xlabel('f/Hz \rightarrow'); ylabel('Yk');
title('3.3 Linienspektrum mit BLACKMAN-Fenster',... 
  'FontWeight','bold'); 

pause;
close all;
figure;
N=512;
k=[0:N-1];
y=100*sin(pi*1980*k/fa)+2*sin(pi*2020*k/fa)+2*rand(1,N)-1;
y=y.*boxcar(N)';
Yk=abs(fft(y));
Ymax=max(Yk); %Maximum feststellen
Yk=20*log10(Yk/Ymax);
f=k*fa/N;
plot(f,Yk,'bo-'); grid; axis([0 fa/2 -100 0]);
xlabel('f/Hz \rightarrow'); ylabel('Yk');
title('3.3 Linienspektrum mit BOXCAR-Fenster',... 
  'FontWeight','bold'); 

pause;
close all;
figure;
N=512;
k=[0:N-1];
y=100*sin(pi*1980*k/fa)+2*sin(pi*2020*k/fa)+2*rand(1,N)-1;
y=y.*chebwin(N)';
Yk=abs(fft(y));
Ymax=max(Yk); %Maximum feststellen
Yk=20*log10(Yk/Ymax);
f=k*fa/N;
plot(f,Yk,'bo-'); grid; axis([0 fa/2 -100 0]);
xlabel('f/Hz \rightarrow'); ylabel('Yk');
title('3.3 Linienspektrum mit CHEBWIN-Fenster',... 
  'FontWeight','bold'); 

pause;
close all;
figure;
N=512;
k=[0:N-1];
y=100*sin(pi*1980*k/fa)+2*sin(pi*2020*k/fa)+2*rand(1,N)-1;
y=y.*hamming(N)';
Yk=abs(fft(y));
Ymax=max(Yk); %Maximum feststellen
Yk=20*log10(Yk/Ymax);
f=k*fa/N;
plot(f,Yk,'bo-'); grid; axis([0 fa/2 -100 0]);
xlabel('f/Hz \rightarrow'); ylabel('Yk');
title('3.3 Linienspektrum mit HAMMING-Fenster',... 
  'FontWeight','bold'); 

pause;
close all;
figure;
N=512;
k=[0:N-1];
y=100*sin(pi*1980*k/fa)+2*sin(pi*2020*k/fa)+2*rand(1,N)-1;
y=y.*kaiser(N)';
Yk=abs(fft(y));
Ymax=max(Yk); %Maximum feststellen
Yk=20*log10(Yk/Ymax);
f=k*fa/N;
plot(f,Yk,'bo-'); grid; axis([0 fa/2 -100 0]);
xlabel('f/Hz \rightarrow'); ylabel('Yk');
title('3.3 Linienspektrum mit KAISER-Fenster',... 
  'FontWeight','bold'); 

pause;
close all;
figure;
N=512;
k=[0:N-1];
y=100*sin(pi*1980*k/fa)+2*sin(pi*2020*k/fa)+2*rand(1,N)-1;
y=y.*triang(N)';
Yk=abs(fft(y));
Ymax=max(Yk); %Maximum feststellen
Yk=20*log10(Yk/Ymax);
f=k*fa/N;
plot(f,Yk,'bo-'); grid; axis([0 fa/2 -100 0]);
xlabel('f/Hz \rightarrow'); ylabel('Yk');
title('3.3 Linienspektrum mit TRIANG-Fenster',... 
  'FontWeight','bold'); 

pause;
close all;
clear all;
clc;

echo on;
%% --- 3.4 Matlab Animation (Ev. Demo in der Vorlesung):
%*********************************************************************
echo off;

% Movie zur FFT Analyse
n=65; fa=4000; N=128; k=[0:N-1]; t=k/fa; f=k*fa/N;
%
f0=fa/4; %Startfrequenz
for j=1:n
fs=f0 +fa/(N*10)*(j-1);
y=sin(2*pi*fs*t)+0.00001*rand(1,N);
%
% hier Fenster wahlweise auskommentieren
y=y.*hann(N,'periodic')';
%y= y.*hamming(N,'periodic')';
Yk=1/N*abs(fft(y));
Ymax=max(Yk); %Maximum feststellen
Yk=20*log10(Yk);
plot(f,Yk,'bo-'); grid;
axis([0 fa/2 -100 0]);
M(j) = getframe;
end
movie(M,1,5) % movie wird zum Schluss 3 Mal gespielt mit 5 Frames per second

pause;
close all;
clear all;
clc;

echo on;
%*********************************************************************
% 4 Kurzzeit-Spektralanalyse
%*********************************************************************
echo off;

%Beispiel-Code:
fa= 8000; ta=1/fa; N=2^14 ; time=[0:N-1]*ta; % Abtastfrequenz und Zeitraster für N Werte
ft=fa/4; %Trägerfrequenz in der Mitte des Basisbereichs
df=fa/N; f=[0:N-1]*df ; % Frequenzvektor

% Bestimmung einer geeigneten Signalfrequenz fs:
% Da das erwartete Spektrum ein Bessellinienspektrum um ft herum
% mit Linien im Abstand fs sein soll, ist es sinnvoll die vorkommenden
% Frequenzen auf Auswertefrequenzen zu legen, d.h. fs ist ein ganzzahliges
% Vielfaches von df
fs=1; % Vorgabe Signalfrequenz
fs=df*round(fs/df) % Auf Auswertefrequenzen angepasste Signalfrequenz
%
% Signalvektor erzeugen:
etha=400; % Modulationsindex
hub=etha*fs % Daraus resultierender Hub; einseitig!!
yfm=cos(2*pi*ft*time + etha*cos(2*pi*fs*time));
% Betrags-Spektrum
Yspek=abs(fft(yfm)); Yspek=Yspek/max(Yspek);
%
% Zweigeteilten Plot erzeugen
subplot(2,1,1); plot(time,yfm); subplot(2,1,2); plot(f,Yspek);
%
% jetzt Ermittlung der Kurzzeitspektren jeweis über
% eine Segmentlänge von N2=512 Werten;
N2=128; w2=hanning(N2); % Fensterfunktion für Kurzzeitspektren
i2=0:N2-1; df2=fa/N2; f2=i2*df2;
% Zeitlicher Überlapp der verschiedenen Fenster M= N2/2;
M=N2/2;
% somit können
Nseg=N/(N2/2) -1; % Segmente ausgewertet werden
%
Y3D=zeros(Nseg, N2/2); % Speicherplatz reservieren!!
%
% Schleife über Nseg Segmente;
for k=1:Nseg;
i_min=(k-1)*M+1; i_max=i_min+N2-1; % Start - Stop - Index
YY= fft(w2'.*yfm(i_min:i_max));
Y3D(k,:)=abs(YY(1:M)); % nur pos. Basisbereich weiterverarbeiten
end;
time2=[0:Nseg-1]*M*ta;
Ynorm=max(max(Y3D)); % maximum suchen
Y3D=Y3D/Ynorm; f2b=f2(1:M);
figure; h=waterfall(f2b,time2,Y3D); view(30,30); %Wasserfall-Plot
%
figure; specgram(yfm,N2,fa); % Matlab Internes Spectogramm

