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

 

Website News

13.08.2017: New Gas e.V. sections (e.g. pictures or prices) added.

05.08.2017: The first "slide carousel" added.

03.08.2017: Started with the RED project. Will be ready in a few years ...

30.07.2017: List of PhD theses updated and extended.

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

 

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

Jens Reermann Defended his Dissertation with Distinction

On Friday, 21st of June, Jens Reermann defended his research on signals processing for magnetoelectric sensor systems very successfully. After 90 minutes of talk and question time he finished his PhD with distinction. Congratulations, Jens, from the entire DSS team.

Jens worked for about three and a half years - as part of the collaborative research center (SFB) 1261 - on all kinds of signal ...


Read more ...