Quelltexte zur Vorlesung Signalverarbeitung

Signalausgabe

Einfache Routine zur Signalausgabe, falls die Routine sound aus octave auf dem Rechner nicht funktioniert. Benötigt den Befehl play z.B. aus dem Paket SoX.

%% my_sound.m
function []=my_sound(snd,Fs)
if(nargin == 1)
Fs = 8000 % MATLAB's default
end
if(rows(snd) != length(snd))
snd = snd';
end
file = tmpnam();
file =[file, '.wav'];

wavwrite( snd, Fs, file );

system(['play ' file]);
system(['rm ' file]);

Abspielen einer (diskreten) Sinuswelle von 440 Hz (Kammerton a)

% play_sine.m
Fs = 8000;
snd = sin( ( 1:Fs )/Fs * 2 * pi * 440 );
my_sound( snd );

Minimalversion von wvtool

MATLABs wvtool bietet die Möglichkeit, ein System in Frequenz sowie Zeitachse darzustellen. Dieser Code ist eine Minimalversion davon.

%% my_wvtool.m
% input: Impulsantwort des Filters
function myfig=my_wvtool(filter)
figure()
subplot(1,2,1);
plot(filter);
title('Zeitachse');
xlabel('Abtastpunkte');
ylabel('Amplitude');

subplot(1,2,2);
N=1024
freq=abs(fft(filter,N));
logfreq=20*log10(freq(1:512));
plot([1:512]/512,logfreq);grid;
title('Frequenzachse');
ylabel('Magnitude (dB)')
xlabel('normalisierte Frequenz');
end

% test my_wvtool
filter = [ 1 4 6 4 1 ];
my_wvtool( filter );

Berechnung des Betragsfrequenzgangs in einem vorgegebenen Frequenzbereich

filterB = [ 1 4 6 4 1 ]; % FIR
filterA = [ 1 ]; % IIR
N = 1000 % Abtastpunkte im Frequenzbereich
Fn = 8000; % Abtastrate
f_max = 50; % maximale Frequenz
f_min = 100; % minimale Frequenz
f = (0:N)/N * (f_max-f_min) + f_min; % Berechnung von N+1 diskreten Frequenzen
e = exp( sqrt(-1) * f/Fn * 2 * pi ); % Berechnung der Punkte auf dem Einheitskreis
H = polyval( filterB, e ) ./ polyval( filterA, e ); % Auswertung von H(e^j omega)
plot( abs( H ) ); % Darstellung des Betragsfrequenzgangs an den Abtastpunkten

50 Hz Notchfilter

%% notch50.m
Fs = 8000;
Fnyquist = Fs/2;
Fnotch = 50;
Frel = Fnotch/Fnyquist;

Q = 10;
width = Frel/Q;
fzeros = [ exp( -sqrt(-1)*pi*Frel ), exp( sqrt(-1)*pi*Frel ) ];
fpoles = (1-width) * fzeros;

% Singnal mit den Frequenzen 30Hz, 40Hz, 50Hz, 60Hz
snd = sin( ( 1:Fs )/Fs * 2 * pi * 30 ) + ...
sin( ( 1:Fs )/Fs * 2 * pi * 40 ) + ...
sin( ( 1:Fs )/Fs * 2 * pi * 50 ) + ...
sin( ( 1:Fs )/Fs * 2 * pi * 60 );

% Erstelle Filter als Polynom in z aus den Null und Polstellen
b = poly(fzeros);
a = poly(fpoles);

% Wende Filter auf Testsignal an
out = filter(b, a, snd );

% Frequenzen in der Eingabe
figure();
freqz(snd);

% Frequenzen in der Ausgabe
figure();
freqz( out );

LPC

N=10 % Length of Filter
S=41 % Length of Signal

% setup signal
sig = 0:(S-1)
sig2 = sin(sig./10*2*pi) + sin(sig./20*2*pi) + sin(sig./40*2*pi)
plot(sig2)

% setup system of equations
M=zeros(S,S)
for i=1:(S-1); M(i,:)=circshift( fliplr(sig2), [0,i]); end

% the current value
b=M(:,1);

% the 5 previous values
M=M(:,2:2+N-1)

% solve Ma = b
a=inverse((transpose(M)*M)) * transpose(M) * b

sol=zeros(1,S);
% test solution by backwards multiplication
for i=1:S;
tmp=circshift(sig2,[0,i-1]);
sol(1,S+1-i) = tmp(1:N) * a;
end

% plot result
plot(transpose([sig2;sol]))

% roots(a)

Web Site Revisions

2013-06-24: initial release