Lecture "Signals and Systems II"

Central limit theorem

%**************************************************************************
% 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

Matlabdemo for exercise ''Zweiseitenbandmodulation''

%--------------------------------------------------------------------------
% Signals and Systems II - Exercise 1
%
% Sound example for the "encoded" speech signal.
%
%
% Nov. 2012, Nov. 2013, Nov. 2014
% Digital Signal Processing and System Theory
% Faculty of Engineering Christian-Albrechts-Universität zu Kiel
% Jochen Withopf, 
%--------------------------------------------------------------------------

clear

%--------------------------------------------------------------------------
% Set parameters, initialize variables
%--------------------------------------------------------------------------
fs = 100000; % has to be high enough for the modulated signals

% --- Minimum and maximum signal frequency
fmin =  300;
fmax = 3400;

% --- Modulation frequencies
f1    = 20000;
f2    = 25000;

% --- Filter cutoff frequencies (3dB cutoffs, we don't use stop band and
%     pass band edge frequencies for simplicity here) and filter order N
fc_hp = 20000;
fc_tp = 10000;
N     = 12;

% --- Design butterworth/Chebychev filters
% [b_hp a_hp] = butter( N, fc_hp/fs*2, 'high');
% [b_tp a_tp] = butter( N, fc_tp/fs*2, 'low');
[b_hp, a_hp] = cheby2( N, 30, fc_hp/fs*2, 'high');
[b_tp, a_tp] = cheby2( N, 30, fc_tp/fs*2, 'low');

% --- Load input signal, bandlimit and upsample it to sample rate fs;
%[u, fs_sig] = wavread( 'signals/interview_heino.wav' );
[u, fs_sig] = wavread( 'signals/sprachnotiz_1.wav' );

[b, a] = butter( 4, [fmin fmax]/fs_sig*2 );
u = filter( b, a, u );
u = resample( u, fs, fs_sig );

%--------------------------------------------------------------------------
% Modulation system
%--------------------------------------------------------------------------
t = linspace(0, length(u)/fs, length(u));
t = t(:);    % make t a column vector

g1 = u.*sin(2*pi*f1*t);
g2 = filter(b_hp, a_hp, g1);
g3 = g2.*sin(2*pi*f2*t);
g4 = filter(b_tp, a_tp, g3);

%--------------------------------------------------------------------------
% Plots and sound
%--------------------------------------------------------------------------

N_win = 1024*4;

[pu, f] = pwelch(u, N_win,128,N_win,fs);
pg1     = pwelch(g1,N_win,128,N_win,fs);
pg2     = pwelch(g2,N_win,128,N_win,fs);
pg3     = pwelch(g3,N_win,128,N_win,fs);
pg4     = pwelch(g4,N_win,128,N_win,fs);

% [pu, f] = pwelch(u,N_win,128,N_win,fs,'twosided');
% pg1    = pwelch(g1,N_win,128,N_win,fs,'twosided');
% pg2    = pwelch(g2,N_win,128,N_win,fs,'twosided');
% pg3    = pwelch(g3,N_win,128,N_win,fs,'twosided');
% pg4    = pwelch(g4,N_win,128,N_win,fs,'twosided');

figure(1);
plot(f/1000, 10*log10([pu pg1 pg2 pg3 pg4]), 'LineWidth', 1);
legend('|U(f)|', '|G_1(f)|', '|G_2(f)|', '|G_3(f)|', '|G_4(f)|');
grid on;
xlabel('Frequency in kHz');
ylabel('Magnitude in dB');

%--------------------------------------------------------------------------
% Spectrograms
%--------------------------------------------------------------------------
u_dn  = resample(u, fs_sig, fs);
g4_dn = resample(g4, fs_sig, fs);
M     = round(40/1000*fs_sig);
ovl   = round(M*0.9);
t_dn  = linspace(0, length(u_dn)/fs_sig, length(u_dn));

figure(2);
sp(1) = subplot(4,2,1);
plot(t_dn, u_dn);
title('Eingangssignal u(t)');
grid on;
xlim([t_dn(1) t_dn(end)]);

sp(2) = subplot(4,2,2);
plot(t_dn, g4_dn);
title('Ausgangssignal g_4(t)');
grid on;

sp(3) = subplot(4,2,[3,5,7]);
[S,F,T] = spectrogram(u_dn, M, ovl, M, fs_sig);
imagesc(T,F,20*log10(abs(S)));
axis xy
ylim([0 fs_sig/4]);
xlabel('Zeit in Sekunden');
ylabel('Frequenz in Hz');

sp(4) = subplot(4,2,[4, 6, 8]);
[S,F,T] = spectrogram(g4_dn, M, ovl, M, fs_sig);
imagesc(T,F,20*log10(abs(S)));
axis xy
ylim([0 fs_sig/4]);
xlabel('Zeit in Sekunden');
ylabel('Frequenz in Hz');

linkaxes(sp, 'x');

Matlabdemo for exercise ''Amplitudenmodulation''

%% 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 as a zip file

  Matlab examples from the exercises