%**************************************************************************
% Matlab introduction
% SuS II (WS21)
%**************************************************************************

close all              % Close all plots
clear                   % Clear workspace
clc                      % Clear console

% Set all fonts to latex style
set(groot, 'DefaultTextInterpreter', 'LaTeX');
set(groot, 'DefaultAxesTickLabelInterpreter', 'LaTeX');
set(groot, 'DefaultAxesFontName', 'LaTeX');
set(groot, 'DefaultAxesFontSize', 14);
set(groot, 'DefaultLegendInterpreter', 'LaTeX'); 

% Save color map in variable
cmap = lines;

% Get info from matlab docu e.g for command 'close' -------------
% help close % Short overview
% doc close  % Detailed documentation

% Matrix and element wise arithmetic operations ---------------------
% Use ; to suppress output of statement
% a = [1, 2, 3, 4];    % Row vector
% a = [1 2 3 4];       % Same row vector
% b = [2; 3; 4; 5];    % Column vector
% a = a';                 % Transpose vector
% c = b * a;             % Matrix multiplication
% c = a.* a;             % Element wise multiplication
% size(c)                 % Get matrix dimensions
% length(c)             % Get max of matrix dim.


%**************************************************************************
% Basic parameters
%**************************************************************************
f_s = 100;            % Sample rate [Hz]
T = 10;                 % Signal duration [s]
f_0 = 0.5;             % Test signal frequency [Hz]

%**************************************************************************
% Auxiliary parameters
%**************************************************************************
L = f_s * T + 1;     % Signal length 

% Generate equidistant time scale from 0 to 10 s
t = (0 : L - 1) / f_s;

%**************************************************************************
% Generate signals
%**************************************************************************
% Generate sine test signal with frequency f_0
v = sin(2*pi*f_0*t); 

% Generate uniformly distributed noise signal
v_u = rand(1, L);

% Generate gaussian noise signal
v_g = randn(1, L);

%**************************************************************************
% Create plots
%**************************************************************************

%%
% Minimal plot of test signal ----------------------------------------------
figure            % Open new figure
plot(t, v_g)        % Actual plotting operation

%%
% Plot with some beautification--------------------------------------------
figure                          
plot(t, v, 'LineWidth', 2, 'Color', 'red')           
title('Test signal with some beautification')
xlabel('Time / s')          % X axis label
ylabel('Magnitude')      % Y axis label
legend('Test signal')    % Add legend
grid on                        % Add legend
xlim([0 1/f_0])              % Limit x-axis to one period

%%
% Discrete stem plot ------------------------------------------------------
figure                          
stem(t, v)           
title('Test signal with stem plot')
xlabel('Time / s')          % X axis label
ylabel('Magnitude')      % Y axis label
legend('Test signal')    % Add legend
grid on                        % Add legend
xlim([0 0.5])                 % Limit x-axis


%%
% Multiple plots with hold on ----------------------------------------------
figure
plot(t, v)
title('Multiple plots with hold on')
hold on                 % Keep previous plot
plot(t, -v)
xlabel('Time / s')                                  
ylabel('Magnitude')                    
legend('Test signal', 'Inv. test signal') 
grid on                                              
xlim([0 1/f_0])                                   

%%
% Multiple plots with subplot ----------------------------------------------
% 2 vertical,1 horizontal plot
figure
subplot(2, 1, 1)    % First subplot
plot(t, v)      
title('Multiple plots with subplot')
ylabel('Magnitude')                    
legend('Test signal') 
grid on 
xlim([0 1/f_0])    

subplot(2, 1, 2)    % Second subplot
plot(t, -v)
xlabel('Time / s')                                  
ylabel('Magnitude')                    
legend('Inv. test signal') 
grid on                                              
xlim([0 1/f_0])    

%**************************************************************************
% 3D Plots (will be discussed next week)
%**************************************************************************

%%
% Plot v, v_2 and v_3 as 3D trajectory ----------------------------------

% Generate sine signal as first component
v_x = sin(2*pi*f_0*t).*t; 
% Generate cosine signal as second component
v_y = cos(2*pi*f_0*t).*t; 
% Generate linear function as third component
v_z = t;

figure                         
plot3(v_x, v_y, v_z)              % 3D plot
title('3D plot as trajectory')
xlabel('X / m')   
ylabel('Y / m')       
zlabel('Z / m')  
axis image                           % Proportional axis
grid on
%%
% Plot noise as surface plot ----------------------------------------------

% Generate gaussian noise as 2D matrix
L_2 = 50;                              % Use less samples
v_2d = randn(L_2, L_2);        % Generate noise
x =  1 : L_2;                           % X axis
y =  1 : L_2;                           % Y axis

% Surface plot
figure      
subplot(2, 1, 1)
surf(x, y, v_2d, 'EdgeColor', 'none', 'FaceColor', 'interp')
title('Surface plot 3D')
c = colorbar;                     % Use colorbar               
ylabel(c, 'Magnitude')
colormap jet
xlabel('X / m')   
ylabel('Y / m')       
zlabel('Z / m')  

subplot(2, 1, 2)
surf(x, y, v_2d, 'EdgeColor', 'none', 'FaceColor', 'interp')
title('Surface plot topview')
c = colorbar;
ylabel(c, 'Magnitude')
xlabel('X / m')   
ylabel('Y / m')       
zlabel('Z / m')  
view(0, 90)                        % Rotate view to XY
axis tight                           

%**************************************************************************
% Amplitude modulation and signal obfuscation
% SuS II (WS21)
%**************************************************************************

close all hidden
clear
clc

set(groot, 'DefaultTextInterpreter', 'LaTeX');
set(groot, 'DefaultAxesTickLabelInterpreter', 'LaTeX');
set(groot, 'DefaultAxesFontName', 'LaTeX');
set(groot, 'DefaultAxesFontSize', 14);
set(groot, 'DefaultLegendInterpreter', 'LaTeX'); 

%**************************************************************************
% Basic parameters
%**************************************************************************
f_s      = 120e3;        % Sampling frequency in Hz
f_low    = 2e3;        % Lowest signal frequency in Hz
f_high   = 4e3;        % Highest signal frequency in Hz
f_mod1   = 20e3;         % Modulation frequency in Hz
f_mod2   = 25e3;         % Modulation frequency 2 in Hz
T_signal = 10;           % Signal length in s

%**************************************************************************
% Calculate support variables
%**************************************************************************
bw = f_high - f_low;                           % Bandwidth in Hz
L_signal = T_signal * f_s;                     % Signal length in samples
t = (0 : L_signal - 1) / f_s;                  % Time scale
f = linspace(-f_s/2, f_s/2, L_signal)./1e3;    % Frequency scale
v_mod1 = sin(2*pi*f_mod1*t)';                  % Modulation signal
v_mod2 = sin(2*pi*f_mod2*t)';                  % Modulation signal 2

%**************************************************************************
% Design filters
%**************************************************************************
% Design FIR highpass filter (complicated...)

f_hp = f_mod1 + [-1e3, 1e3];
a_hp = [0 1]; 
rp = 1;
rs = 120;

dev = [(10^(rp/20)-1)/(10^(rp/20)+1)  10^(-rs/20)]; 
[n_ord_bp,fo_bp,ao_bp,w] = firpmord(f_hp,a_hp,dev,f_s);

b_high = firpm(n_ord_bp+1,fo_bp,ao_bp);


% Design Butterworth lowpass filter for demodulation -----------------------------
f_c = 30e3;      % Cutoff freq in Hz
n = 30;          % Cutoff freq in Hz

[z,p,k] = butter(n, f_c/(f_s/2));
[sos_low, g_low] = zp2sos(z, p, k);

%**************************************************************************
% Generate initial time signal
%**************************************************************************
% Use sum of multiple sine functions in specified frequency range

u = 1* sin(2*pi*f_low*t)';
u = u + 0.5 *sin(2*pi*f_high*t)';
%%
%**************************************************************************
% Processing chain for amplitude modulation
%**************************************************************************
% Modulate, demodulate and filter inital signal
% Use either this or the other processing chains

% Process FFT of initial time signal --------------------------------------
% Process FFT:                  fft(..)
% Absolute spectrum:            abs(..)
% Weight correctly:             ../ L
% Shift to get correct x axis:  circshift(.., - ceil(L / 2)
% See MATLAB docu for more info on FFT processing 
U = circshift(abs(fft(u)) / L_signal, - ceil(L_signal/2));

% Modulation step ---------------------------------------------------------
g_1 = u.*v_mod1;
G_1 = circshift(abs(fft(g_1)) / L_signal, - ceil(L_signal/2));

% Demodulation step -------------------------------------------------------
g_2 = g_1.*v_mod1;
G_2 = circshift(abs(fft(g_2)) / L_signal, - ceil(L_signal/2));

% Lowpass filter step -----------------------------------------------------
g_3 = sosfilt(sos_low, g_2);
G_3 = circshift(abs(fft(g_3)) / L_signal, - ceil(L_signal/2));

% Not used ----------------------------------------------------------------
g_4 = zeros(L_signal, 1);
G_4 = zeros(L_signal, 1);

%%
%**************************************************************************
% Processing chain for signal obfuscation 
%**************************************************************************
% Modulate, highpass filter, modulate again, lowpass filter
% Insert "u = g_4;" in second turn to reconstruct initial signal

% Process FFT -------------------------------------------------------------
U = circshift(abs(fft(u)) / L_signal, - ceil(L_signal/2));

% Modulation step ---------------------------------------------------------
g_1 = u.*v_mod1;
G_1 = circshift(abs(fft(g_1)) / L_signal, - ceil(L_signal/2));

% Supress lower frequency band --------------------------------------------
%g_2 = sosfilt(sosHigh, g_1);
g_2 = filtfilt(b_high, 1, g_1);
G_2 = circshift(abs(fft(g_2)) / L_signal, - ceil(L_signal/2));

% Modulation step ---------------------------------------------------------
g_3 = g_2.*v_mod2;
G_3 = circshift(abs(fft(g_3)) / L_signal, - ceil(L_signal/2));

% Lowpass filter step -----------------------------------------------------
g_4 = filtfilt(sos_low, g_low, g_3); %filter(b_high, 1, g_3)
G_4 = circshift(abs(fft(g_4)) / L_signal, - ceil(L_signal/2));

%**************************************************************************
% Processing chain for signal obfuscation reversal
%**************************************************************************
% Modulate, highpass filter, modulate again, lowpass filter

% Modulation step ---------------------------------------------------------
g_5 = g_4.*v_mod1;
G_5 = circshift(abs(fft(g_5)) / L_signal, - ceil(L_signal/2));

% Supress lower frequency band --------------------------------------------
g_6 = filtfilt(b_high, 1, g_5);
G_6 = circshift(abs(fft(g_6)) / L_signal, - ceil(L_signal/2));

% Modulation step ---------------------------------------------------------
g_7 = g_6.*v_mod2;
G_7 = circshift(abs(fft(g_7)) / L_signal, - ceil(L_signal/2));

% Lowpass filter step -----------------------------------------------------
g_8 = filtfilt(sos_low, g_low, g_7); %filter(b_high, 1, g_3)
G_8 = circshift(abs(fft(g_8)) / L_signal, - ceil(L_signal/2));

%%
%**************************************************************************
% Plots
%**************************************************************************
% Evaluate section by section to get plots after each other

t_range = 1e-2; % Limited time axis scale is easier to view

figure
set(gcf, 'Position', [2000 0 1400 1400]) % Set window size

subplot(5, 2, 1)
plot(t, u)
xlim([0, t_range])
ylabel('\(u(t)\)')
grid on

subplot(5, 2, 2)
plot(f, U)
ylabel('\(|U_(jw)|\)')
grid on

%%
subplot(5, 2, 3)
plot(t, g_1)
xlim([0, t_range])
ylabel('\(g_\mathrm{1}(t)\)')
grid on

subplot(5, 2, 4)
plot(f, G_1)
ylabel('\(|G_\mathrm{1}(jw)|\)')
grid on

%%
subplot(5, 2, 5)
plot(t, g_2)
xlim([0, t_range])
ylabel('\(g_\mathrm{2}(t)\)')
grid on

subplot(5, 2, 6)
plot(f, G_2)
ylabel('\(|G_\mathrm{2}(jw)|\)')
grid on

%%
subplot(5, 2, 7)
plot(t, g_3)
xlim([0, t_range])
ylabel('\(g_\mathrm{3}(t)\)')
grid on

subplot(5, 2, 8)
plot(f, G_3)
ylabel('\(|G_\mathrm{3}(jw)|\)')
grid on

%%
subplot(5, 2, 9)
plot(t, g_4)
xlim([0, t_range])
ylabel('\(g_\mathrm{4}(t)\)')
xlabel('\(t\) / s')
grid on

subplot(5, 2, 10)
plot(f, G_4)
ylabel('\(|G_\mathrm{4}(jw)|\)')
xlabel('\(\frac{\omega}{2\pi}\) / kHz')
grid on

%%
%**************************************************************************
% Plots
%**************************************************************************
% Evaluate section by section to get plots after each other

t_range = 1e-2; % Limited time axis scale is easier to view

figure
set(gcf, 'Position', [2000 0 1400 1400]) % Set window size

subplot(5, 2, 1)
plot(t, u)
xlim([0, t_range])
ylabel('\(u(t)\)')
grid on

subplot(5, 2, 2)
plot(f, U)
ylabel('\(|U_(jw)|\)')
grid on

%%
subplot(5, 2, 3)
plot(t, g_5)
xlim([0, t_range])
ylabel('\(g_\mathrm{5}(t)\)')
grid on

subplot(5, 2, 4)
plot(f, G_5)
ylabel('\(|G_\mathrm{5}(jw)|\)')
grid on

%%
subplot(5, 2, 5)
plot(t, g_6)
xlim([0, t_range])
ylabel('\(g_\mathrm{6}(t)\)')
grid on

subplot(5, 2, 6)
plot(f, G_6)
ylabel('\(|G_\mathrm{6}(jw)|\)')
grid on

%%
subplot(5, 2, 7)
plot(t, g_7)
xlim([0, t_range])
ylabel('\(g_\mathrm{7}(t)\)')
grid on

subplot(5, 2, 8)
plot(f, G_7)
ylabel('\(|G_\mathrm{7}(jw)|\)')
grid on

%%
subplot(5, 2, 9)
plot(t, g_8)
xlim([0, t_range])
ylabel('\(g_\mathrm{8}(t)\)')
xlabel('\(t\) / s')
grid on

subplot(5, 2, 10)
plot(f, G_8)
ylabel('\(|G_\mathrm{8}(jw)|\)')
xlabel('\(\frac{\omega}{2\pi}\) / kHz')
grid on

%**************************************************************************
% Angle modulation
% SuS II (WS21)
%**************************************************************************

close all hidden
clear
clc

cmap = lines;

set(groot, 'DefaultTextInterpreter', 'LaTeX');
set(groot, 'DefaultAxesTickLabelInterpreter', 'LaTeX');
set(groot, 'DefaultAxesFontName', 'LaTeX');
set(groot, 'DefaultAxesFontSize', 14);
set(groot, 'DefaultLegendInterpreter', 'LaTeX'); 

%**************************************************************************
% Basic parameters
%**************************************************************************
f_s         = 500;          % Sampling frequency in Hz
f_mod    = 20;            % Modulation frequency in Hz
f_0         = 1;             % Frequency of initial signal in Hz
T_signal = 10;           % Signal length in s

%**************************************************************************
% Calculate support variables
%**************************************************************************
L_signal = T_signal * f_s;                     % Signal length in samples
w_0 = 2*pi*f_0;                                    
w_mod = 2*pi*f_mod;
f = linspace(-f_s/2, f_s/2, L_signal);    % Frequency scale
t = ((0 : L_signal - 1) / f_s)';                % Time scale

%**************************************************************************
% Processing chain for different analog modulation schemes
%**************************************************************************
% Initial time signal and spectrum -----------------------------
v = cos(2*pi*f_0*t);                          
V = circshift(abs(fft(v)) / L_signal, - ceil(L_signal/2));

% Amplitude modulated signal and spectrum --------------
v_am = v.*cos(w_mod*t);
V_am = circshift(abs(fft(v_am)) / L_signal, - ceil(L_signal/2));

% Frequency modulation -----------------------------------------
% Instantaneous frequency and phase
k = w_0;        % Use k to adapt frequency shift                                    
Omega_fm = w_mod + k*2*pi*v;
Phi_fm = w_mod*t + k*2*pi*cumsum(v)*1/f_s;

% Modulated signal and spectrum
v_fm = cos(Phi_fm);
V_fm = circshift(abs(fft(v_fm)) / L_signal, - ceil(L_signal/2));

% Phase modulation -----------------------------------------------
% Instantaneous frequency and phase
k = 1;
Omega_pm = w_mod + k*2*pi*diff([v; v(end)])*f_s;
Phi_pm = w_mod*t + k*2*pi*v;

% Modulated signal and spectrum
v_pm = cos(Phi_pm);
V_pm = circshift(abs(fft(v_pm)) / L_signal, - ceil(L_signal/2));

%**************************************************************************
% Plots
%**************************************************************************
% Evaluate section by section to get plots after each other

tRange = [0, 2 * 1/f_0];   % Limited time axis scale is easier to view
fRange = 3 * [-f_mod, f_mod]; 

figure
%set(gcf, 'Position', [2000 0 1400 1400]) % Set window size

% Initial signal -----------------------------
subplot(6, 2, 1)
plot(t, v)
xlim(tRange)
ylabel('\(v(t)\)', 'FontSize', 20)
xlabel('\(t\) [s]', 'FontSize', 20);
title('Input - Time signal', 'FontSize', 18)
grid on

%%
subplot(6, 2, 2)
plot(f, V)
xlim(fRange)
ylabel('\(|V_(jw)|\)', 'FontSize', 20)
xlabel('\(\frac{\omega}{2\pi}\)  [Hz]', 'FontSize', 20);
title('Input - Spectrum', 'FontSize', 18)
grid on

%%
% Amplitude modulation -----------------------------
subplot(6, 2, 3)
plot(t, v_am, 'Color', cmap(2, :))
xlim(tRange)
ylabel('\(v_\mathrm{am}(t)\)', 'FontSize', 20)
xlabel('\(t\) [s]', 'FontSize', 20);
title('DSB - Time signal', 'FontSize', 18)
grid on

%%
subplot(6, 2, 4)
plot(f, V_am, 'Color', cmap(2, :))
xlim(fRange)
ylabel('\(|V_\mathrm{am}(jw)|\)', 'FontSize', 20)
xlabel('\(\frac{\omega}{2\pi}\)  [Hz]', 'FontSize', 20);
title('DSB - Spectrum', 'FontSize', 18)
grid on

%%
% Frequency modulation -----------------------------
subplot(6, 2, 5)
plot(t, Omega_fm, 'Color', cmap(3, :))
xlim(tRange)
ylabel('\(\Omega_\mathrm{T,fm}(t)\) / \(\mathrm{s^{-1}}\)', 'FontSize', 20)
xlabel('\(t\) [s]', 'FontSize', 20);
title('Frequency modulation - Instan. frequency', 'FontSize', 18)
grid on

%%
subplot(6, 2, 6)
plot(t, Phi_fm, 'Color', cmap(3, :))
xlim(tRange)
ylabel('\(\Phi_\mathrm{T,fm}(t)\) / rad', 'FontSize', 20)
xlabel('\(t\) [s]', 'FontSize', 20);
title('Frequency modulation - Instan. phase', 'FontSize', 18)
grid on

%%
subplot(6, 2, 7)
plot(t, v_fm, 'Color', cmap(3, :))
xlim(tRange)
ylabel('\(v_\mathrm{fm}(t)\)', 'FontSize', 20)
xlabel('\(t\) [s]', 'FontSize', 20);
title('Frequency modulation - Time signal', 'FontSize', 18)
grid on

%%
subplot(6, 2, 8)
plot(f, V_fm, 'Color', cmap(3, :))
xlim(fRange)
ylabel('\(|V_\mathrm{fm}(jw)|\)', 'FontSize', 20)
xlabel('\(\frac{\omega}{2\pi}\)  [Hz]', 'FontSize', 20);
title('Frequency modulation - Spectrum', 'FontSize', 18)
grid on

%%
% Phase modulation -----------------------------
subplot(6, 2, 9)
plot(t, Omega_pm, 'Color', cmap(4, :))
xlim(tRange)
ylabel('\(\Omega_\mathrm{T,pm}(t)\) / \(\mathrm{s^{-1}}\)', 'FontSize', 20)
xlabel('\(t\) [s]', 'FontSize', 20);
title('Phase modulation - Instan. frequency', 'FontSize', 18)
grid on

%%
subplot(6, 2, 10)
plot(t, Phi_pm, 'Color', cmap(4, :))
xlim(tRange)
ylabel('\(\Phi_\mathrm{T,pm}(t)\)  / rad', 'FontSize', 20)
xlabel('\(t\) [s]', 'FontSize', 20);
title('Phase modulation - Instan. phase', 'FontSize', 18)
grid on

%%
subplot(6, 2, 11)
plot(t, v_pm, 'Color', cmap(4, :))
xlim(tRange)
ylabel('\(v_\mathrm{pm}(t)\)', 'FontSize', 20)
xlabel('\(t\) [s]', 'FontSize', 20);
title('Phase modulation - Time signal', 'FontSize', 18)
grid on

%%
subplot(6, 2, 12)
plot(f, V_pm, 'Color', cmap(4, :))
xlim(fRange)
ylabel('\(|V_\mathrm{pm}(jw)|\)', 'FontSize', 20)
xlabel('\(\frac{\omega}{2\pi}\)  [Hz]', 'FontSize', 20);
title('Phase modulation - Spectrum', 'FontSize', 18)
grid on

%**************************************************************************
% Pulse amplitude modulation
% SuS II (WS21)
%**************************************************************************

close all hidden
clear
clc

cmap = lines;

set(groot, 'DefaultTextInterpreter', 'LaTeX');
set(groot, 'DefaultAxesTickLabelInterpreter', 'LaTeX');
set(groot, 'DefaultAxesFontName', 'LaTeX');
set(groot, 'DefaultAxesFontSize', 14);
set(groot, 'DefaultLegendInterpreter', 'LaTeX'); 

%**************************************************************************
% Basic parameters
%**************************************************************************
f_s         = 500;          % Sampling frequency in Hz
f_mod    = 20;            % Modulation frequency in Hz
f_0         = 1;             % Frequency of initial signal in Hz
T_signal = 10;           % Signal length in s
f_steep = 5;             % Transition band width in Hz

%**************************************************************************
% Calculate support variables
%**************************************************************************
L_signal = T_signal * f_s;                     % Signal length in samples
T_mod = 1 / f_mod;
f = linspace(-f_s/2, f_s/2, L_signal);    % Frequency scale
t = ((0 : L_signal - 1) / f_s)';                % Time scale

%**************************************************************************
% Processing chain for pulse amplitude modulation
%**************************************************************************
% Initial time signal and spectrum -----------------------------
v = cos(2*pi*f_0*t);                          
V = circshift(abs(fft(v)) / L_signal, - ceil(L_signal/2));

% Impulse train carrier ---------------------------------------------
v_mod = zeros(length(t), 1);
v_mod(~mod(t, T_mod)) = v_mod(~mod(t, T_mod)) + 1;
V_mod = circshift(abs(fft(v_mod)) / L_signal, - ceil(L_signal/2));

v_mod2 = cos(2*pi*f_mod*t);

% Amplitude modulated signal and spectrum --------------
v_pam = v.*v_mod;
V_pam = circshift(abs(fft(v_pam)) / L_signal, - ceil(L_signal/2));

v_am = v.*v_mod2;
V_am = circshift(abs(fft(v_am)) / L_signal, - ceil(L_signal/2));
V_am = V_am / max(V_am) * max(V_pam);        % Adjust magnitude

% Simple butterworth lowpass filter for demodulation ------
f_c = 10;      % Cutoff freq in Hz
n = 10;         % Filter order
[z,p,k] = butter(n, f_c/(f_s/2));
sosLow = zp2sos(z, p, k);

% FIR lowpass filter for (improved) demodulation ------------
f_lp = f_mod/2 + [-1 1] * f_steep;
a_lp = [1 0]; 
rp = 1;
rs = 120;

dev = [(10^(rp/20)-1)/(10^(rp/20)+1)  10^(-rs/20)]; 
[n_ord_lp,fo_lp,ao_lp,w] = firpmord(f_lp,a_lp,dev,f_s);
b_lp = firpm(n_ord_lp+1,fo_lp,ao_lp);

% Plot frequency response if desired
% freqz(b_lp, 1, 2^16, f_s);

% Lowpass filter step ----------------------------------------------------------------
% Choose Butterworth filter...
% v_pam_lp = sosfilt(sosLow, v_pam);
% [h, f_2] =  freqz(sosLow, 2^10 * 10, f_s);

% ...or FIR filter
v_pam_lp = filter(b_lp, 1, v_pam);
[h, f_2] =  freqz(b_lp, 1, 2^10 * 10, f_s);


% Compute filter frequency response
f_2 = [-flip(f_2); f_2];
h = abs(h);
h = [flip(h); h];
h = h * max(V_pam); % Adjust magnitude

V_pam_lp = circshift(abs(fft(v_pam_lp)) / L_signal, - ceil(L_signal/2));

%**************************************************************************
% Plots
%**************************************************************************
% Evaluate section by section to get plots after each other

tRange = [0, 3 * 1/f_0];   % Limited time axis scale is easier to view
fRange = 3 * [-f_mod, f_mod]; 

figure
%set(gcf, 'Position', [2000 0 1400 1400]) % Set window size

% Initial signal -----------------------------
subplot(4, 2, 1)
plot(t, v)
xlim(tRange)
ylabel('\(v(t)\)', 'FontSize', 20)
xlabel('\(t\) [s]', 'FontSize', 20);
grid on
title('Input - Time signal', 'FontSize', 18)

%%
subplot(4, 2, 2)
plot(f, V)
xlim(fRange)
ylabel('\(|V_(jw)|\)', 'FontSize', 20)
xlabel('\(\frac{\omega}{2\pi}\) [Hz]', 'FontSize', 20);
grid on
title('Input - Spectrum', 'FontSize', 18)

%%
% Impulse train carrier-----------------
subplot(4, 2, 3)
plot(t, v_mod, 'Color', cmap(1, :))
xlim(tRange)
ylabel('\(v_\mathrm{c}(t)\)', 'FontSize', 20)
xlabel('\(t\) [s]', 'FontSize', 20);
grid on
title('Impulse train - Time signal', 'FontSize', 18)

%%
subplot(4, 2, 4)
plot(f, V_mod, 'Color', cmap(1, :))
xlim(fRange)
ylabel('\(|V_\mathrm{c}(jw)|\)', 'FontSize', 20)
xlabel('\(\frac{\omega}{2\pi}\) [Hz]', 'FontSize', 20);
grid on
title('Impulse train - Spectrum', 'FontSize', 18)

%%
% PAM modulated signal-----------------
subplot(4, 2, 5)
plot(t, v_pam, 'Color', cmap(1, :))
xlim(tRange)
ylabel('\(v_\mathrm{pam}(t)\)', 'FontSize', 20)
xlabel('\(t\) [s]', 'FontSize', 20);
grid on
hold on
title('Modulation - Time signal', 'FontSize', 18)

%%
subplot(4, 2, 6)
plot(f, V_pam, 'Color', cmap(1, :))
xlim(fRange)
ylabel('\(|V_\mathrm{pam}(jw)|\)', 'FontSize', 20)
xlabel('\(\frac{\omega}{2\pi}\) [Hz]', 'FontSize', 20);
grid on
hold on
title('Modulation - Spectrum', 'FontSize', 18)

%%
% Amplitude modulated signal-----------------
subplot(4, 2, 5)
plot(t, v_am, 'LineStyle', '-', 'LineWidth', 0.2)
legend('PAM', 'AM')

subplot(4, 2, 6)
plot(f, V_am, 'LineStyle', '-', 'LineWidth', 0.2)

%%
% Lowpass filter frequency response
subplot(4, 2, 6)
plot(f_2, h, 'Color', cmap(3, :), 'LineStyle', '--', 'LineWidth', 2)
legend('PAM', 'AM', 'Lowpass')

%%
% Lowpass demodulation -----------------------------
subplot(4, 2, 7)
plot(t, v_pam_lp, 'Color', cmap(1, :))
xlim(tRange)
ylabel('\(v_\mathrm{lp}(t)\)', 'FontSize', 20)
xlabel('\(t\) [s]', 'FontSize', 20);
grid on
title('Demodulation - Time signal', 'FontSize', 18)

%%
subplot(4, 2, 8)
plot(f, V_pam_lp, 'Color', cmap(1, :))
xlim(fRange)
ylabel('\(|V_\mathrm{lp}(jw)|\)', 'FontSize', 20)
xlabel('\(\frac{\omega}{2\pi}\) [Hz]', 'FontSize', 20);
grid on
title('Demodulation - Spectrum', 'FontSize', 18)

%**************************************************************************
% Parameter fuer die Generierung von gleichverteilten Zufallsvariablen
%**************************************************************************
m_v   =      0.0;        % Mittelwert
var_v =      1.0;        % Varianz
N     = 400000;          % Anzahl der Variablen pro Experiment
M     =      4;          % Anzahl der summierten Zufallsvariablen

%**************************************************************************
% Erzeugen der Zufallszahlen
%**************************************************************************
V = ((rand(M,N) - 0.5) / sqrt(12) * sqrt(var_v)) + m_v;

%**************************************************************************
% Nicht-normierte Histogramm-Analyse der ersten Zufallsvariablen
%**************************************************************************
fig = figure(1);
set(fig,'Units','Normalized');
set(fig,'Position',[0.1 0.1 0.8 0.8]);

[h,x] = hist(V(1,:),20);
bar(x,h,0.7);
ylabel('Haeufigkeit');
xlabel('Variable v_1');
grid on

%return;

%**************************************************************************
% Histogramm-Analyse der einzelnen Zufallsvariablen
%**************************************************************************
A_min    =   -0.6;
A_max    =    0.6;
IntWidth =    (A_max - A_min) / 50;

IntVec   = A_min:IntWidth:A_max;
HistRes  = zeros(M,length(IntVec));

for k = 1:M
    HistRes(k,:) = hist(V(k,:),IntVec) / (N * IntWidth);
end;

%**************************************************************************
% Darstellung der einzelnen Histogrammanalysen
%**************************************************************************
fig = figure(2);
set(fig,'Units','Normalized');
set(fig,'Position',[0.1 0.1 0.8 0.8]);

for k = 1:min(4,M)
   subplot(2,2,k)
   bar(IntVec,HistRes(k,:),0.7)
   grid on;
   xlabel(['Variable v_',num2str(k)']);
   axis([A_min A_max -0.05*max(HistRes(k,:)) 1.05*max(HistRes(k,:))]);
   ylabel('Relative Haeufigkeit');   
end;

%return;

%**************************************************************************
% Summation der Zufallsvariablen
%**************************************************************************
v = zeros(1,N);
for k = 1:M
    v = v + V(k,:);
end;

%**************************************************************************
% Histogrammanalyse der Summation
%**************************************************************************
histRes = hist(v,IntVec) / (N * IntWidth);

%**************************************************************************
% Darstellung der Histogrammanalyse der Summe
%**************************************************************************
fig = figure(3);
set(fig,'Units','Normalized');
set(fig,'Position',[0.1 0.1 0.8 0.8]);

bar(IntVec,histRes,0.7)
grid on;
xlabel('Variable y');
axis([A_min A_max -0.05*max(histRes) 1.05*max(histRes)]);
ylabel('Relative Haeufigkeit');   

%**************************************************************************
% Darstellung einer Gauss-Dichte
%**************************************************************************
f_g = 1/(std(v)*sqrt(2*pi)) * exp( -0.5 * ((IntVec-mean(v))/std(v)).^2);
hold on
plot(IntVec,f_g,'r','LineWidth',2);
hold off

%% Example for single-sideband modulation and "double"-sideband modulation
%
% -> general concept
% -> advantage of singe-sibeband modulation (reduced error caused by frequency shift)
%
% Jens Reermann
%% Reset
clear all;
close all;
clc;
%% Parameter
f_s=44.1e3;
PLAY_SIGNAL=0;
AUDIO_ENABLED = 0;

%%
fc=16e3;
f_error=1500;
phi_error = 0; %pi*0.48;

%% Load signal
if AUDIO_ENABLED == 1
    [u, fs_sig] = audioread( 'signals/interview_heino.wav' );
    n=0:length(u)-1;
else
    fsig=2000;
    fs_sig=100e3;
    T=10;
    Nsig=T*fs_sig;
    n=0:Nsig-1;
    u = sin(2*pi*fsig/fs_sig.*n)'+randn(Nsig,1);
end


%% LP 3kHz | 4kHz
f_l=3e3;
f_h=4e3;
f = [ f_l f_h];
a = [1 0];
rp = 1;           
rs = 120;     
dev = [ (10^(rp/20)-1)/(10^(rp/20)+1) 10^(-rs/20)]; 
[n_ord,fo,ao,w] = firpmord(f,a,dev,f_s);
b = firpm(n_ord+1,fo,ao);
u_f = filter(b,1,u);

%% Carrier signals

c=cos(2*pi*n*fc/f_s)';
c_e=cos(2*pi*n*(fc+f_error)/f_s+phi_error)';

%% Modulation

% DSB
m=c.*u_f;

% BP: DSB->SSB
f_bp = [ fc fc+1000 fc+f_l fc+f_h];
a_bp = [0 1  0]; 
dev_bp = [10^(-rs/20) (10^(rp/20)-1)/(10^(rp/20)+1) 10^(-rs/20)]; 
[n_ord_bp,fo_bp,ao_bp,w] = firpmord(f_bp,a_bp,dev_bp,f_s);
b_bp = firpm(n_ord_bp+1,fo_bp,ao_bp);

m_bp = filtfilt(b_bp,1,m);    %SSB

%% Demodulation

% *cos
u_d=m.*c;
u_d_bp=m_bp.*c;
u_de=m.*c_e;
u_de_bp=m_bp.*c_e;
% LP
u_df = filter(b,1,u_d);
u_df_bp = filter(b,1,u_d_bp);
u_def = filter(b,1,u_de);
u_def_bp = filter(b,1,u_de_bp);

%% Play signalss
if PLAY_SIGNAL==1
    sound(u_f,f_s)
    pause
    sound(u_df,f_s)
    pause
    sound(u_def,f_s)
    pause
    sound(u_def_bp,f_s)
end

%% PSD estimation
N_FFT=2^12;
[Puu, F] = pwelch(u,hann(N_FFT),N_FFT-N_FFT/4,N_FFT,f_s);
[Pufuf, F] = pwelch(u_f,hann(N_FFT),N_FFT-N_FFT/4,N_FFT,f_s);
[Pmm, F] = pwelch(m,hann(N_FFT),N_FFT-N_FFT/4,N_FFT,f_s);
[Pmm_bp, F] = pwelch(m_bp,hann(N_FFT),N_FFT-N_FFT/4,N_FFT,f_s);
[Pudfudf, F] = pwelch(u_df,hann(N_FFT),N_FFT-N_FFT/4,N_FFT,f_s);
[Pudefudef, F] = pwelch(u_def,hann(N_FFT),N_FFT-N_FFT/4,N_FFT,f_s);
[Pudfudf_bp, F] = pwelch(u_df_bp,hann(N_FFT),N_FFT-N_FFT/4,N_FFT,f_s);
[Pudefudef_bp, F] = pwelch(u_def_bp,hann(N_FFT),N_FFT-N_FFT/4,N_FFT,f_s);

%% Plotting (Modulation)
[val1 ind1] =max(10*log10(Puu));
[val2 ind2] =max(10*log10(Pmm));
[val3 ind3] =max(10*log10(Pmm_bp));

if AUDIO_ENABLED
    y_lims=[-160 -50];
else
    y_lims=[-100 10];
end

figure

plot(F, 10*log10(Puu));
text(F(ind1),val1,[num2str(val1) 'dB']);
legend('S_{uu}');
xlabel('f / Hz');ylabel('LDS / dB');xlim([0 f_s/2]);ylim(y_lims);
pause
hold on;
plot(F, 10*log10(Pufuf),'r');
legend('S_{uu}','S_{ufuf}');
xlabel('f / Hz');ylabel('LDS / dB');xlim([0 f_s/2]);ylim(y_lims);
pause
plot(F, 10*log10(Pmm),'g');
text(F(ind2),val2,[num2str(val2) 'dB']);
legend('S_{uu}','S_{ufuf}','S_{mm}');
xlabel('f / Hz');ylabel('LDS / dB');xlim([0 f_s/2]);ylim(y_lims);
pause
plot(F, 10*log10(Pmm_bp),'k');
text(F(ind3),val3,[num2str(val3) 'dB']);
legend('S_{uu}','S_{ufuf}','S_{mm}','S_{mm_{BP}}');
xlabel('f / Hz');ylabel('LDS / dB');xlim([0 f_s/2]);ylim(y_lims);
pause
hold off;

%% Plotting (Demodulation)

[val3 ind3] =max(10*log10(Pufuf));
[val4 ind4] =max(10*log10(Pudfudf));
[val5 ind5] =max(10*log10(Pudefudef));
[val6 ind6] =max(10*log10(Pudefudef_bp));

figure
xlim([0 f_h+1000]);
ylabel('LDS / dB');
xlabel('f / Hz');
if AUDIO_ENABLED
    ylim=[-160 -50];
else
    ylim=[-100 10];
end

hold on;
plot(F, 10*log10(Pufuf));
text(F(ind3),val3,[num2str(val3) 'dB']);
legend('S_{ufuf}');
pause
plot(F, 10*log10(Pudfudf),'g');
text(F(ind4),val4,[num2str(val4) 'dB']);
legend('S_{ufuf}','S_{ZSB}');
pause
plot(F, 10*log10(Pudefudef),'r');
text(F(ind5),val5,[num2str(val5) 'dB']);
legend('S_{ufuf}','S_{ZSB}','S_{eeZSB}');
pause
plot(F, 10*log10(Pudefudef_bp),'k');
text(F(ind6),val6,[num2str(val6) 'dB']);
legend('S_{ufuf}','S_{ZSB}','S_{eeZSB}','S_{eeESB}');
pause

hold off;

disp( [ num2str(20*log10(cos(phi_error))) ] );

  Matlab examples from the exercises