iirnotch in matlab

[Sebastian H. (sebihain)]

hallo forum,

ich muss ein digitales notchfilter (50hz, fs=500hz) erstellen und 
berechnen. Habe mich hier im Forum umgesehen und einige sehr gute sachen 
gefunden. daraus  hab ich mir mal folgendes gebaut (meine 
matlabkentnisse sind sehr beschränkt!):

function test
fs=1000; % doppelte von F-Abtast
%
t = 0:.001:.239; % beliebiger Zeitraum
x = sin(2*pi*50*t)+sin(2*pi*150*t)+sin(2*pi*200*t);
%
[vFrequency,vAmplitude]=test_fft(x,fs);  % FFT mit Hamming Window

w0=50/250;
bw=w0/35;
[num,den]=iirnotch(w0,bw);
x_n=filter(num,den,x);
[f_notch,amp_notch]=test_fft(x_n,fs);

subplot(4,1,1);plot(t,x);title('gestörtes Signal')
subplot(4,1,2);plot(vFrequency, vAmplitude);title('FFT - gestörtes 
Signal')
subplot(4,1,3);plot(t,x_n);title('notch-filter Signal')
subplot(4,1,4);plot(f_notch,amp_notch);title('FFT - notch Signal')
return


Aber leider bekomme ich die 50hz nicht weg!!! Kann mir jemand sagen was 
ich da falsch gemacht habe???

hoffe auf eure hilfe (und hoffe ich hab nichts vergessen zu erwähnen)

vielen Dank
Sebi

[Detlef _. (detlef_a)]

Du nimmst 1kHz samplerate, obwohl Du in der Einleitung von 500Hz 
sprichst. w0 ist auch (vermutlich) bezogen auf 500Hz. Die Nullstellen 
von num/den ausrechnen, dann weiß Du Bescheid!

Cheers
Detlef

[Sebastian H. (sebihain)]

Hallo Detlef,

ja da mit der Abtastf. ist mir ein fehler unterlaufen, aber wenn ich es 
ändere macht das filter trotzdem nichts. ich bekomme für num und den je 
drei werte. die setze ich dann brav wie matlab es sagt in 
(x_n=filter(num,den,x)) ein, aber die 50hz f. taucht immer noch auf.

w0 = f0/(f_abtast/2) und muss zwischen 0 und 1 liegen
bw = gibt die bandbreite des notch an <1

oder ist das falsch???  ich versteh gar nix mehr.

hoffe du kannst mir weiterhelfen.
Danke Sebi

[Detlef _. (detlef_a)]

Die Nullstellen von num und den ausrechnen: roots(num). Da kommen 
jeweils zwei konjugiert komplexe Werte raus. Für richtige num und den 
ist der Winkel dieser beiden Werte +/-(50/1000)*pi. Der Betrag von num 
muß eins sein, der Betrag von den muß gegen 1 gehen für kleiner werdende 
Bandbreite.

Statt iirnotch einfach so:
num = poly([exp(j*50*pi/1000) exp(-j*50*pi/1000)]);
betrag =0.97;
den = poly([betrag*exp(j*50*pi/1000) betrag*exp(-j*50*pi/1000)]);

Cheers
Detlef

[Tommi H. (drmota)]

Das Notchfilter  funktioniert. Mit

t = (0:1/1024:0.1-1/1024);          %Abtastzeit 1/1024 s
x = sin(2*pi*50*t)+sin(2*pi*150*t)+sin(2*pi*200*t);
y = fft(x);
m = abs(y);
f = (0:length(y)-1)*1023/length(y);
plot(f,y);
w0=50/512;
bw=w0/5;            %ist mit 5 schneller als mit 35
[num,den]=iirnotch(w0,bw);


num,den habe ich dann in ein Simulinkmodell verfrachtet. Dieses  Modell
zeigt (siehe Anhang) dass nach etwa 0.4s die Filterung stattfindet.

[Tommi H. (drmota)]

Hier noch das Simulinkmodell.

[Sebastian H. (sebihain)]

Vielen Dank für eure Antworten,

hat jetzt geklappt. hab aber noch ein verständnisproblem.

w0 wird ja wie folgt definiert --> w0 = f0/(f_abtast/2)

dabei ist f0 die Frequenz die ich nicht haben will und f_abtast eben die 
abtastfrequenz. wenn ich jetzt aber mit w0=50/(500/2) das ganze laufen 
lasse funktioniert es nicht, aber wenn ich mit w0=50/500 setze 
funktioniert es???

und bw gibt mit die Bandweite an? je kleiner bw wird, desto 
wahrscheinlicher ist es, dass er die frequenz nicht mehr ganz filtern 
kann.

Sind diese Überlegungen richtig?? oder alles schachsinn?? kann mir da 
vielleicht jemand licht in mein dunkles hirn bringen?

vielen dank
Sebi

PS: kann jemand ein buch empfehlen das sich mit DSP beschäfftigt - und 
es von anfang an erklärt.(im bestfall noch dazu einfach!!!)   -   DANKE