Hallo, der kontinuierliche Fall ist klar: y(t)=x(t)***g(t) mit g(t) Impulsantwort des Systems, x(t) Eingangssignal und *** als Faltungsoperation. ausgeschrieben mit Integrall: y(t)=integral von -inf bis +inf über g(tau)*x(t-tau)*dtau ist das System kausal, dann gilt: g(t<0)=0 das es gilt kann man sich so veranschaulichen, dass die Gewichtung nur von x(t) Werten vorliegt, die bereits in der Vergangenheit liegen einschliesslich des Jetzt. Also kann y(t) vereinfacht werden zu: y(t)=integral von 0 bis +inf über g(tau)*x(t-tau)*dtau Wenn x(t) auch für negative Werte 0 ist, dann lässt sich es nochmal vereinfachen: y(t)=integral von 0 bis t über g(tau)*x(t-tau)*dtau Darauf ansetzend die Formel für den diskreten Fall: y[n]=Summe i=0 bis n über g[i]*x[n-i] g[i] und x[n-i] müssen aber denselben Zeitpunkt repräsentieren. Wie sieht es denn jetzt konkret im diskreten Fall aus? Angenommen ich nehme einfachen passiven RC Filter mit g(t)=1/(R*C)*exp(-t/(R*C)) Jetzt habe ich ein diskretes x[n] Signal, mit einer Samplerate Ts abgetastet. Ich würde also denken, dass ich g(t) ebenfalls mit Ts samplen muss. Ist das korrekt? Wenn ja, dann würde ja die Wahl von Ts nicht nur von der höchsten im x(t) vorhandenen Frequenz abhängen, sondern auch von der höchsten vorhandenen Frequenz in g(t). fsample=2*max(max(fx),max(fg)) also. Danke im Voraus.
:
Verschoben durch Admin
mein Versuch in Octave und Check in LTspice zeigt eine unverkennbare Übereinstimmung.
1 | t=0:0.001:1; |
2 | N=size(t)(2); |
3 | x=sin(2*pi*5*t)+sin(2*pi*15*t)+sin(2*pi*100*t); |
4 | |
5 | function y=faltung(g,x,k) |
6 | tmp=0; |
7 | for i=1:k; |
8 | tmp+=g(i)*x(k-i+1); # DSP MAC - Operation? |
9 | end; |
10 | y=tmp; |
11 | end; |
12 | |
13 | R=100; |
14 | C=0.001; |
15 | RC=R*C; |
16 | g=1/RC*exp(-t/RC); |
17 | |
18 | y=[]; |
19 | for k=1:N; |
20 | y(k)=faltung(g,x,k); |
21 | end; |
22 | |
23 | figure(1); |
24 | clf; |
25 | subplot(3,1,1); |
26 | plot(t,x); |
27 | grid on; |
28 | xlabel 'time'; |
29 | ylabel 'input'; |
30 | subplot(3,1,2); |
31 | plot(t,g); |
32 | grid on; |
33 | xlabel 'time'; |
34 | ylabel 'g impulse response'; |
35 | subplot(3,1,3); |
36 | plot(t,y); |
37 | grid on; |
38 | xlabel 'time'; |
39 | ylabel 'signal response'; |
40 | |
41 | pause; |
1) Ist "tmp+=g(i)*x(k-i+1)" Zeile die berühmt berüchtigte MAC - Operation? 2) Impulseantwort ist infinite, dh wird erst im Unendlichen 0. Dennoch habe ich das Gefühl, dass meine Softimplementierung eine FIR ist. Diese verwendet keine y Rückkopplung zur Berechnung von nachkommenden y Werten. Aber dadurch, dass viele Koeffizienten verwendet werden, ist das Ergebnis nicht schlecht. Wie kann ich IIR Variante implementieren? Grüsse, Daniel
Hallo! Daniel -------- schrieb: > Wie sieht es denn jetzt konkret im diskreten Fall aus? Angenommen > ich nehme einfachen passiven RC Filter mit g(t)=1/(R*C)*exp(-t/(R*C)) > Jetzt habe ich ein diskretes x[n] Signal, mit einer Samplerate Ts > abgetastet. Ich würde also denken, dass ich g(t) ebenfalls mit Ts > samplen muss. Ist das korrekt? Warum willst Du denn eine andere Abtastfrequenz nehmen? Und warum willst Du ausgerechnet die Impulsantwort eines analogen TP abtasten? Wenn Du eh schon mit zeitdiskreten Werten arbeitest gibt es doch wesentlich besseres um zu Filtern. Oder soll das nur ein Beispiel sein? Daniel -------- schrieb: > 1) Ist "tmp+=g(i)*x(k-i+1)" Zeile die berühmt berüchtigte MAC - > Operation? Ja ;) > 2) Impulseantwort ist infinite, dh wird erst im Unendlichen 0. > Dennoch habe ich das Gefühl, dass meine Softimplementierung eine FIR > ist. Ja das ist ein FIR-Filter. Der FIR-Filter macht ja auch nichts anderes, als das Eingangssignal mit seiner Impulsantwort zu falten. Grüsse Loup
zuerst kleine Korrektur an dem Octave Skript: global N=1000;
1 | function y=faltung(g,x,k) |
2 | global N; |
3 | tmp=0; |
4 | for i=1:k; |
5 | tmp+=g(i)*x(k-i+1); |
6 | end; |
7 | y=tmp/N; # <-- divide by N |
8 | end; |
ansonsten stimmen die Werte nicht mit der simulierten Version überein. @Loup >Warum willst Du denn eine andere Abtastfrequenz nehmen? >Und warum willst Du ausgerechnet die Impulsantwort eines analogen TP >abtasten? Wenn Du eh schon mit zeitdiskreten Werten arbeitest gibt es >doch wesentlich besseres um zu Filtern. Oder soll das nur ein Beispiel >sein? Am Anfang habe ich ja noch keine zeitdiskreten Werte für g[i]. Ich habe nur die g(t) Funktion. Der zeitliche Abstand zwischen g[i],g[i+1],... Koeffizienten bedeutet doch nichts anderes als Abtastung. Es ist aus meiner Sicht nicht die Frage warum ich abtasten muss, sondern wie. Und da denke ich, dass hier ebenso das Nyquist Kriterium greift, wie bei ganz normalen Signalabtastung. Ich habe keinen strengen Beweis dazu, dass das so sein muss. Wenn nun g(t) sich aus höheren Frequenzen zusammensetzt, dann bekomme ich g[i],g[i+1] mit einem kleineren zeitlichen Abstand, als für x[i],x[i+1]. Dann kann ich nicht g[i]*x[k-i] einfach so multiplizieren, weil die Zeitskalierungen nicht identisch sind. Also müsste x(t) auch so abgetastet werden wie g(t). Darum wollte ich zuerst die maximale Frequenz aus den beiden Funktionen bestimmen. >Ja das ist ein FIR-Filter. Der FIR-Filter macht ja auch nichts anderes, >als das Eingangssignal mit seiner Impulsantwort zu falten. Ok, dann habe ich richtig vermutet :) Weisst Du wie ich IIR Filter für diesen "einfachen" Fall mit g(t)=-1/RC*exp(-t/RC) bauen kann? Ich weiss nicht genau, aber glaube, dass wenn g(t) erst im Unendlichen verschwindet, so wie bei uns, dann sollte auch eine IIR Variante aufstellbar sein. Grüsse
Daniel -------- schrieb: > Und da denke ich, dass hier ebenso das Nyquist Kriterium greift, wie > bei ganz normalen Signalabtastung. Ich habe keinen strengen Beweis dazu, > dass das so sein muss. Wenn die Abtastfrequenz die selbe ist hast Du jedenfalls am wenigsten Ärger. Ohne das jetzt genau durchdacht zu haben würde ich sagen, dass es sonst ein paar seltsame Effekte geben könnte, wenn Du ein Signal mit Aliasing mit einem anderen faltest. Bei Deiner speziellen Impulsantwort mit der abfallenden e-Funktion geht es aber auch nach der Nyquistfreqenz erst mal relativ konstant (klein) weiter bis die (gespiegelte) e-Funktion wieder bis zur Abtastfrequenz hin ansteigt. Sie ist quasi ihr eigenes Anti-Aliasing-Filter. Ändere doch mal Dein Programm entsprechend ab und probiere aus was passiert. Würde mich auch interessieren! Daniel -------- schrieb: > Weisst Du wie ich IIR Filter für diesen "einfachen" Fall mit > g(t)=-1/RC*exp(-t/RC) bauen kann? Ich weiss nicht genau, aber > glaube, dass wenn g(t) erst im Unendlichen verschwindet, so wie > bei uns, dann sollte auch eine IIR Variante aufstellbar sein. Ja kann man schon machen, aber es ist dann keine Faltung mehr. Korrespondenz für die z-Transformation:
Führt auf die Differenzengleichung:
Der Betrag von a muss natürlich kleiner 1 sein sonst wirds instabil. Siehe dazu auch: http://www.dspguide.com/ch19/2.htm Grüsse
ich bin jetzt etwas weiter ;-) Mir ist im Nachhinein klar geworden, dass ich mich wahrscheinlich falsch ausgedrückt habe. Wenn g(t) im Unendlichen 0 wird, dann muss die Implementierung es zwangsläufig approximieren, wenn man den Filter als Summe von Gewichtungen macht. Nachteilig ist, dass man viele Koeffizienten brauchen wird => viel Platz in FPGA(etc) Oder man implementiert es rekursiv => platzsparender. Wahrscheinlich (im Allgmeinen kann ich es nicht sagen) wird es sich dabei auch um eine Approximaion handeln. Ich habe hier einen weiteren Entwurf, diesmal rekursiv.
1 | % octave/matlab script |
2 | |
3 | % using the model |
4 | % |
5 | % 1 |
6 | % G(s) = ------ |
7 | % RC*s+1 |
8 | |
9 | % using bilinear approximation |
10 | % |
11 | % 2 1-z^{-1} |
12 | % s = - * -------- |
13 | % T 1+z^{-1} |
14 | |
15 | % after the substitution is carried out, we get |
16 | % |
17 | % y(i)*(RC*2/T+1) = x(i)+x(i-1)-y(i-1)*(1-RC*2/T); |
18 | % |
19 | |
20 | R=10; |
21 | C=0.1; |
22 | RC=R*C; |
23 | T=0.01; |
24 | x=[zeros(1,1000),ones(1,1000)]; |
25 | N=size(x)(2); |
26 | y=[0]; |
27 | for i=2:2000; |
28 | y(i)=(1/(RC*2/T+1))*(x(i)+x(i-1)-y(i-1)*(1-RC*2/T)); |
29 | end; |
30 | |
31 | subplot(3,1,1); |
32 | plot(linspace(0,T*N,N),x); |
33 | grid; |
34 | |
35 | subplot(3,1,2); |
36 | plot(linspace(0,T*N,N),y); |
37 | grid; |
38 | set(gca,'xtick',0:1:T*N); |
39 | |
40 | % verifing the result |
41 | m=tf(1,[RC 1]); |
42 | subplot(3,1,3); |
43 | plot(linspace(0,T*N,N), lsim(m,x',linspace(0,T*N,N))); |
44 | grid; |
45 | set(gca,'xtick',0:1:T*N); |
46 | |
47 | pause; |
wegen der Vollständigkeit auch für den Hochpass
1 | % octave/matlab script |
2 | |
3 | % using the model |
4 | % |
5 | % RC*s |
6 | % G(s) = ------ |
7 | % RC*s+1 |
8 | |
9 | % using bilinear approximation |
10 | % |
11 | % 2 1-z^{-1} |
12 | % s = - * -------- |
13 | % T 1+z^{-1} |
14 | |
15 | % after the substitution is carried out, we get |
16 | % |
17 | % y(i)*(RC*2/T+1) = RC*2/T(x(i)+x(i-1))-y(i-1)*(1-RC*2/T); |
18 | % |
19 | |
20 | R=100; |
21 | C=0.01; |
22 | RC=R*C; |
23 | T=0.01; |
24 | x=[zeros(1,1000),ones(1,1000)]; |
25 | y=[0]; |
26 | N=size(x)(2); |
27 | for i=2:N; |
28 | y(i)=(RC*2/T*(x(i)-x(i-1))-(1-RC*2/T)*y(i-1))/(RC*2/T+1); |
29 | end; |
30 | |
31 | t=linspace(0,T*N,N); |
32 | |
33 | subplot(3,1,1); |
34 | plot(t,x,'g'); |
35 | set(gca,'xtick',0:1:N*T); |
36 | xlabel 'time'; |
37 | ylabel 'input'; |
38 | grid; |
39 | |
40 | subplot(3,1,2); |
41 | plot(t,y); |
42 | set(gca,'xtick',0:1:N*T); |
43 | xlabel 'time'; |
44 | ylabel 'output'; |
45 | grid; |
46 | |
47 | m=tf([RC 0],[RC 1]); |
48 | subplot(3,1,3); |
49 | plot(t,lsim(m,x',t),'r'); |
50 | set(gca,'xtick',0:1:N*T); |
51 | xlabel 'time'; |
52 | ylabel 'verify'; |
53 | grid; |
54 | |
55 | pause; |
Daniel -------- schrieb: > Mir ist im Nachhinein klar geworden, dass ich mich wahrscheinlich > falsch ausgedrückt habe. Wenn g(t) im Unendlichen 0 wird, dann muss > die Implementierung es zwangsläufig approximieren, wenn man den Filter > als Summe von Gewichtungen macht. Wenn die Impulsantwort nicht endlich ist, dauert die Faltung unter Umständen halt etwas länger ;-) Irgendwann gehen die Werte der e-Funktion aber eh im Quantisierungsrauschen oder der Rundung unter... > Nachteilig ist, dass man viele > Koeffizienten brauchen wird => viel Platz in FPGA(etc) Du musst ja nicht jedem Rechenschritt einen eigenen Hardware-mac spendieren, da gibt es extra Rechenwerke dazu. Und gerade mit einem FPGA sollte die Geschwindigkeit nicht das Problem sein. > Oder man implementiert es rekursiv => platzsparender. Dafür bekommt man aber andere Probleme wie evtl. Instabilität (kann auch gerne mal durch Rundungsfehler passieren), kein linearer Phasengang usw. > Wahrscheinlich > (im Allgmeinen kann ich es nicht sagen) wird es sich dabei auch > um eine Approximaion handeln. Also speziell bei Deiner Funktion denke ich das es keine Approximaion ist, sondern (bis auf die Rechenungenauigkeit) die exakte Abbildung des RC-Gliedes. Du hast es doch oben bereits ausgerechnet!? Aber ehrlich gesagt weis ich immer noch nicht was das ganze eigentlich werden soll!? Wenn Du schon einen digitalen Tiefpass bauen willst, dann kopiere doch nicht einfach ein so popeliges analoges Filter! In Deinem ersten Script mit der Faltung sind das doch 1000 Schleifendurchläufe, oder? Mit der Größenordnung könntest Du auch einen FIR-TP mit mindestens 100 dB (!) Sperrdämpfung haben können und das bei einer Abfallzeit von nur ein paar Hz! Und 1000 Koeffizienten sind schon sehr viel. Ein normaler TP hat eher so Ordnung 50 bis 200. Grüsse
Bitte melde dich an um einen Beitrag zu schreiben. Anmeldung ist kostenlos und dauert nur eine Minute.
Bestehender Account
Schon ein Account bei Google/GoogleMail? Keine Anmeldung erforderlich!
Mit Google-Account einloggen
Mit Google-Account einloggen
Noch kein Account? Hier anmelden.