Code Examples (Signals and Systems - Part 2)

Lecture

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

 

Exercise

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-Universitt 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))) ] );

 

Recent Publications

P. Durdaut, J. Reermann, S. Zabel, Ch. Kirchhof, E. Quandt, F. Faupel, G. Schmidt, R. Knöchel, and M. Höft: Modeling and Analysis of Noise Sources for Thin-Film Magnetoelectric Sensors Based on the Delta-E Effect, IEEE Transactions on Instrumentation and Measurement, published online, 2017

P. Durdaut, S. Salzer, J. Reermann, V. Röbisch, J. McCord, D. Meyners, E. Quandt, G. Schmidt, R. Knöchel, and M. Höft: Improved Magnetic Frequency Conversion Approach for Magnetoelectric Sensors, IEEE Sensors Letters, published online, 2017

 

Website News

18.06.2017: Page about KiRAT news added (also visible in KiRAT).

31.05.2017: Some pictures added.

23.04.2017: Time line for the lecture "Adaptive Filters" added.

13.04.2017: List of PhD theses added.

Contact

Prof. Dr.-Ing. Gerhard Schmidt

E-Mail: gus@tf.uni-kiel.de

Christian-Albrechts-Universität zu Kiel
Faculty of Engineering
Institute for Electrical Engineering and Information Engineering
Digital Signal Processing and System Theory

Kaiserstr. 2
24143 Kiel, Germany

Recent News

Alexej Namenas - A New Guy in the Team

In June Alexej Namenas started in the DSS Team. He will work on real-time tracking algorithms for SONAR applications. Alexej has done both theses (Bachelor and Master) with us. The Bachelor thesis in audio processing (beamforming) and the Master thesis in the medical field (real-time electro- and magnetocardiography). In addition, he has intership erperience in SONAR processing.

We are pretty ...


Read more ...