www.mikrocontroller.net

Forum: Digitale Signalverarbeitung / DSP RekursivI( IIR )Filter mit Matlab entwerfen


Autor: Stefan Feldmann (wasserman2002)
Datum:
Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Hallo an alle,
es geht um RIAA Entzerrung. [sehen Sie das Bild]
Ich will einen IIR Rekursiv Filter 2-Ordnung (Hochpass) mit Matlab 
entwerfen .
Die Koeffizienten der digitalen Übertragungfunktion b(z^-1)/a(z^-1) in
der absteigenden Potenzen von z,sind:

44,1 KHz % Abtastfrequenz
B =[0,02675918611906 -0,04592084787595 0,019212292972391]
A =[1,00000000000000 -0,73845850035973 -0,17951755477430 ]
Fehler +/-0.25dB

Es ist in Form
Out = b0*in[0] + b1*in[-1] + b2*[-2] -a1*out[-1] -a2*out[-2]
wo in [0, -1, -2] sind die aktuellen Eingabe und der vorherigen 2 ;und
aus [-1,-2] sind in den letzten zwei Ausgänge.

http://www.thel-audioworld.de/module/phono/RIAA.htm
Der Ausgang oder Die lösung sollte so aussehen wie das Bild hier
in Pdf.     http://www.hagtech.com/pdf/riaa.pdf

Autor: Christoph Kessler (db1uq) (christoph_kessler)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Die Lösung im PDF bezieht sich aber auf das inverse RIAA-Filter.

Die passive analoge Schaltung aus zwei C und zwei R ist hier gezeigt:
http://www.supertube.de/RIAA/RIAA.htm
aus der Übertragungsgleichung der analogen Schaltung H(s) kann man 
irgendwie auf die digitale H(z) umrechnen, möglicherweise bietet Matlab 
das als Funktion?

Autor: Unit* (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Was ist die Frage?

Autor: Christoph Kessler (db1uq) (christoph_kessler)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Wieso eigentlich Hochpass? Das ist doch eindeutig eine Tiefpass-Kurve 
mit drei Eckfrequenzen. Oder soll es doch ein inverses RIAA-Filter sein?

Autor: Alex (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Umrechnung zwischen z und s:

c2d bzw. d2c

Autor: Stefan Feldmann (wasserman2002)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hallo,

>Was ist die Frage?
wie kann ich einen IIR Rekursiv Filter 2-Ordnung (Hochpass) mit Matlab
entwerfen?.

Autor: Alex (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Warum nimmst du nicht die integrierte Filter Design Toolbox? Dort kannst 
du dir dein Filter in einem GUI zusammenklicken.

Autor: Unit* (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Wenn du die Filterkoeffizienten hast, dann hast du dein Filter. Was 
willst du noch entwerfen???
Oder du willst es realisieren? Oder wat?
Was wäre, wenn du einfach vernünftig beschreiben würdest, was du genau 
willst, anstatt hier 20 threads zu öffnen, und das Forum mit deiner 
Frage vollzustopfen?

> >Was ist die Frage?
> wie kann ich einen IIR Rekursiv Filter 2-Ordnung (Hochpass) mit Matlab
> entwerfen?.

Autor: Stefan Feldmann (wasserman2002)
Datum:
Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
sorry aber der Aufgabe wurde mir nicht vernünftig erklärt,aber jetzt ist 
klar geworden.
Also die Aufgabe ist einfach einen Sinus sweep (20-20KHz) mit meinem 
Filter filtern.
Ich hab das Programm geschrieben,es gibt trozdem fehler ,weiss ich nicht 
wo die sind,
Das ergebniss solte genau umgehehrt wie das Bild ganz oben.
die sollte wie im Bild in pdf.[X-Achse=Logf,y-Achse=Dämfung].
Ich hoffe,es ist jetzt vernünftig beschrieben.
danke an alle.


%Sinus Sweep(Eingangssignal)
%--------------------------
clear
fsample=44100; %Abtastfrequenz
tend=10;
t=0:tend;
dph=20*2*pi*1/fsample+20000*2*pi*(0:1/fsample:tend)/(tend*fsample);
ph=cumsum(dph);
signal=sin(ph);
signal=abs(fft(signal));
figure (1)
plot(signal(20:3000),'.-')
soundsc(signal,44100)
% Übertragungsfunktion vom Filter
%-------------------------------
B =[0.02675918611906  -0.04592084787595  0.019212292972391];%num
A =[1.00000000000000  -0.73845850035973 -0.17951755477430 ];%den
roots(B);
roots(A);
H=freqz(B,A,441001);%ich glaube hier ist das problem
figure (2)
plot(abs(H));
%Ausgang (gefilterter Eingangssignal)
%---------------------------
for s =1:10000  %nur Probe
OUT(s) = signal(1,s)*abs(H(s,1));
end;
figure (3)
plot(OUT)
soundsc(OUT,44100)
return

Autor: Unit* (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
> Das ergebniss solte genau umgehehrt wie das Bild ganz oben.
> die sollte wie im Bild in pdf.

Wieso sollte es umgekehrt aussehen? Ich nehme an, dass du mit dem Sweep 
eine breitbandige Anregung erzielen willst. Wenn man ein breitbandiges 
flaches Anregungsspektrum hat, dann sieht das Spektrum des 
Ausgangsignals genauso aus, wie die Frequenzkarakteristik des Filters. 
Wenn du diese flippen willst, musst du das Anregungssignal mit (1-H) 
filtern!

Ich würde vorschlagen, dass du das System statt dem Sweep einfach mit 
einem Dirac-Stoß anregst (z.B. x = [1 zeros(1,100000)]) und damit das 
(Inverse-)Filter anregst (filter(a-b,b,x)).

Autor: Christoph Kessler (db1uq) (christoph_kessler)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
In Elektor gabs so ein inverses RIAA-Filter schon in den Siebzigern 
untrer dem Namen "RezipRIAA".
Schallplatten werden schon seit Jahrzehnten mit einer Höhenanhebung oder 
Tiefenabsenkung gepresst, weil nur so eine halbe Stunde auf eine Seite 
passt. Die großen Ausschläge im Bassbereich würden einen größeren 
Spurabstand erfordern. Im Tonabnehmer-Vorverstärker wird diese 
Höhenanhebung mit der anfanges gezeigten Tiefpasskurve wieder rückgängig 
gemacht.

Um diese Vorverstärker zu testen (oder um selber Platten herzustellen) 
braucht man das inverse Hochpassfilter. In Reihe zum Verstärker muß 
wieder ein gerader Frequenzgang herauskommen.

Autor: wasserman2002 (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hallo Unit*(gast),
>Ich würde vorschlagen, dass du das System statt dem Sweep einfach mit
einem Dirac-Stoß anregst (z.B. x = [1 zeros(1,100000)]) und damit das
(Inverse-)Filter anregst (filter(a-b,b,x)).
kannst du mir bitte das Programm schreiben,
danke.

Autor: Unit* (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
> kannst du mir bitte das Programm schreiben

Steht doch in meinem letzten Beitrag!

Autor: Stefan Feldmann (wasserman2002)
Datum:
Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Hallo,
das Ergebinss sollte so ausehen wie die Kurve in (Figure 5) in diesem 
PDF.
http://www.hagtech.com/pdf/riaa.pdf
Ich hab gemacht wie (unit *Gast) gesagt hat,sieht aber nicht gut aus !!
so sieht das Programm aus:

x = [1 zeros(1,44100)]; %Dirac Stoss als Eingang
B =[0.02675918611906  -0.04592084787595  0.019212292972391];%num
A =[1.00000000000000  -0.73845850035973 -0.17951755477430 ];%den
roots(B);
roots(A);
H=freqz(B,A,44100);  %Fs=44100Hz
figure (1)
plot(abs(H),'b')
y=filter(A-B,B,x); %Ausgang (gefilterter Eingangssignal)
figure (2)
plot(y)


Mit simulation hat aber gut funktioniert! das Programm unten gehört zu
Simmulation.

clear IN;
clear H;
clear OUT;
OUT=simouty.signals.values;
  OUT=fft(OUT);
IN=simoutx.signals.values;
  IN=fft(IN);
  for s=1:22050;
    H(s) =OUT(s)/IN(:,:,s);%Out und In sollen gleiche grosse haben
  end
  figure(1)
plot(abs(H));
y=-62:1:-25; %wie kann ich meine Kurve im Bereich [-20 bis 20]
             %verschieben und X-Achse mit "Hz"???
Hlog=20*log10(abs(H));
figure (2);
plot(Hlog)
hold on;
plot(50,y); %erste Eckpunkt =f1
plot(500,y); %zweite Eckpunkt=f2
plot(2122,y); %dritte Eckpunkt=f3
hold off;

Autor: Unit* (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
> x = [1 zeros(1,44100)]; %Dirac Stoss als Eingang
> B =[0.02675918611906  -0.04592084787595  0.019212292972391];%num
> A =[1.00000000000000  -0.73845850035973 -0.17951755477430 ];%den

Die nächsten beiden Zeilen sind unnötig:

> roots(B);
> roots(A);


> H=freqz(B,A,44100);  %Fs=44100Hz
> figure (1)
> plot(abs(H),'b')

Mensch! H = B/A -> 1-H = (A-B)/A => y = filter(A-B,A,x);

> y=filter(A-B,B,x); %Ausgang (gefilterter Eingangssignal)
> figure (2)

Wieso transformierst du das Signal nicht???

> plot(y)


Von mir aus kannst du auch freqz statt filter benutzen:

H2 = freqz(A-B,A,44100);
figure
plot(abs(H2))

Autor: Josef (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
"aus der Übertragungsgleichung der analogen Schaltung H(s) kann man
irgendwie auf die digitale H(z) umrechnen, möglicherweise bietet Matlab
das als Funktion?"

siehe Z-Transformation

Autor: Stefan Feldmann (wasserman2002)
Datum:
Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Hallo unit*(gast)
wenn ich so "y=filter(A-B,B,x) "hreibe sttat so"y=filter(B,A,x)"bekomme 
ich einen tiefpass.und ich will Hochpass haben.sehe nochmal PDF.
ich hab jetzt das programm fertiggestelt,sieht so aus:

B =[0.02675918611906  -0.04592084787595  0.019212292972391];%num
A =[1.00000000000000  -0.73845850035973 -0.17951755477430 ];%den
H=freqz(B,A,44100);  %Fs=44100Hz
figure (1)
plot(abs(H),'r')
x = [1 zeros(1,44100)]; %Dirac Stoss als Eingang
y=filter(B,A,x); %Ausgang (gefilterter Eingangssignal)
y=fft(y);
ylog=20*log10(abs(y));
figure (2)
plot(ylog)

die frage jetzt ,wie kann ich meine Kurve um 40db nach oben 
verschieben?? mit Simulation habe ich einen verstärker am Ausgagan 
angeschlossen,wie geht das denn mit M-file, was soll ich zu meinem 
programm noch schreiben?
danke.

Autor: Unit* (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Du hast aber folgendes geschrieben:

> Das ergebniss solte genau umgehehrt wie das Bild ganz oben.

Ich dachte du meinst das Bild im pdf. Ok.

+40dB:

y=filter(100*B,A,x);

oder

y = 100*y;

oder

plot(ylog+40)

oder

...

Autor: Stefan Feldmann (wasserman2002)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
vielen dank, ich hab das probiert,hat funktioniert.

Antwort schreiben

Die Angabe einer E-Mail-Adresse ist freiwillig. Wenn Sie automatisch per E-Mail über Antworten auf Ihren Beitrag informiert werden möchten, melden Sie sich bitte an.

Wichtige Regeln - erst lesen, dann posten!

  • Groß- und Kleinschreibung verwenden
  • Längeren Sourcecode nicht im Text einfügen, sondern als Dateianhang

Formatierung (mehr Informationen...)

  • [c]C-Code[/c]
  • [avrasm]AVR-Assembler-Code[/avrasm]
  • [code]Code in anderen Sprachen, ASCII-Zeichnungen[/code]
  • [math]Formel in LaTeX-Syntax[/math]
  • [[Titel]] - Link zu Artikel
  • Verweis auf anderen Beitrag einfügen: Rechtsklick auf Beitragstitel,
    "Adresse kopieren", und in den Text einfügen




Bild automatisch verkleinern, falls nötig
Bitte das JPG-Format nur für Fotos und Scans verwenden!
Zeichnungen und Screenshots im PNG- oder
GIF-Format hochladen. Siehe Bildformate.
Hinweis: der ursprüngliche Beitrag ist mehr als 6 Monate alt.
Bitte hier nur auf die ursprüngliche Frage antworten,
für neue Fragen einen neuen Beitrag erstellen.

Mit dem Abschicken bestätigst du, die Nutzungsbedingungen anzuerkennen.