Forum: Digitale Signalverarbeitung / DSP / Machine Learning Diskrete Faltung auslegen


von Daniel -. (root)


Lesenswert?

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
von Daniel -. (root)


Angehängte Dateien:

Lesenswert?

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

von Loup (Gast)


Lesenswert?

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

von Daniel -. (root)


Lesenswert?

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

von Loup (Gast)


Lesenswert?

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

von Daniel -. (root)


Angehängte Dateien:

Lesenswert?

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;

von Daniel -. (root)


Angehängte Dateien:

Lesenswert?

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;

von Loup (Gast)


Lesenswert?

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
Noch kein Account? Hier anmelden.