Verfasst am: 27.03.2009, 18:51
Titel: Faltung u. inverse F. im Frequenzraum mit conv() u deconv()?
Hallo,
ich habe drei verschiedene Eingangssignale (versch. Sinus Signale) für ein LTI System generiert und diese dann ganauso wie die Sprungantwort des Systems in den Frequenzraum transformiert. Dazu habe ich den Befehl fft() verwendet. Anschließend habe ich die transformierten Signale je separat mit der ebenfalls transf. Sprungantwort mit dem Befehl conv() gefaltet.
Ein Plot der in den Zeitbereich zurücktransf. Daten zeigt mir den gefalteten Signalverlauf. Sieht OK aus.
Nun wollte ich die Faltung wieder umkehren und transformierte diesmal die gefalteten Signale wieder in den Frequenzraum. Mit dem Befehl deconv(gef. Signal, Sprungantwort) und anschließdender Rücktransformation erhalte für alle meine 3 Signale einen konstanten Wert im Zeitbereich. Wo sind meine ursprünglichen Sinus Signale hin? Wieso funktioniert die inverse Faltung über diesen Weg nicht?
bin auch gerade am Überlegen wieso das so ist.......
Wie sieht denn dein plot aus nachdem du deconv() benutzt hast, also vor ifft???
Sieht der plot ähnlich aus wie der nach der fft() am Anfang?
Vielleicht liegt der Fehler auch darin, dass du NACH der fft() die Convolution benutzt..........
Ich kenne die Prozedur (bisher) nur so, dass man eine fft() erst nach einer Faltung benutzt........
Werde mal weiter darüber nachdenken, wenn mir was einfällt gebe ich dir bescheid.
Ciao
DeusRa
Gast
Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
Verfasst am: 28.03.2009, 14:09
Titel:
Hallo DeusRa,
zunächst einmal vielen Dank für deine Bemühungen. Vermutlich liegt das Problem tatsächlich schon darin, dass ich erst transformiere und danach falte. Leider kann ich das gerade nicht ausprobieren, da ich die Daten nicht parat habe. Werde ich aber im Laufe des WE's testen und dir bescheid geben.
>Wie sieht denn dein plot aus nachdem du deconv() benutzt hast, also vor >ifft??? Sieht der plot ähnlich aus wie der nach der fft() am Anfang?
Wenn das so wäre, könnte man das Ergebnis wenigstens nachvollziehen. Schaue ich mir auch mal an.
Zwei kleine Fragen hätte ich noch: Wäre die Faltung im Frequenzraum (also nach fft) nicht auch über y = x .* h möglich. D. h. über elementweise Multiplikation beider Matrizen? Wenn ja, muss die Sprungantwort zwangsläufig die gleiche Zeitdauer haben wie das Eingangssignal?
Vielen Dank ersteinmal. Ich melde mich vermutlich erst heute Abend wieder, etwa ab 18:00
Zwei kleine Fragen hätte ich noch: Wäre die Faltung im Frequenzraum (also nach fft) nicht auch über y = x .* h möglich. D. h. über elementweise Multiplikation beider Matrizen? Wenn ja, muss die Sprungantwort zwangsläufig die gleiche Zeitdauer haben wie das Eingangssignal?
Also an sich nicht (glaube ich....)......denn eine Faltung keine reine Multiplikation mit den Werten aus deinen Matrizen ist....... an sich ist Convolution ein Intergral über g(x-y)*f(y) (nach dy). Also dieses x-y ist wichtig....also rein Mathematisch. fft() wird wohl so ähnlich sein........wenn ich deine frage richtig verstehe.
Ich melde mich morgen abend auch nochmal, vielleicht ist mir da auch was eingefallen.
Ciao
DeusRa
...
Gast
Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
Verfasst am: 01.04.2009, 16:52
Titel:
Servus DeusRa,
sorry, für die verspätete Antwort. Ich habe mir am WE einen fiesen Magen-Darm-Virus eingefangen. Heute bin ich aber wieder vor Ort gewesen und konnte deinen Hinweis "Faltung mit conv() ausschließlich im Zeitbereich" berücksichtigen. - Nun funktioniert alles einwandfrei.
Danke für die Hilfe
Mark
Gast
Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
Verfasst am: 17.04.2009, 10:52
Titel:
Hallo,
aufgrund der zukünftig hohen Datenmengen möchte ich die Faltung nun doch im Frequenzraum vornehmen. Die elementweise Multiplikation des transformierten Eingangssignals X mit der transformierten Impulsantwort führt nicht zu korrekten Ergebnissen: Y = X .* H
dabei lautet der diskrete Faltungssatz doch genau so:
{y_k} = {(x*h)_k} <--> {H_j}={X_j H_j}
"_k" und "_j" bedeuten die normalerweise tiefgestellten Indices (Latex geht hier nicht, oder?). Links ist die Faltung im Zeitbereich durch den * Operator und rechts das Pendant im Frequenzraum ausgedrückt. Wieso ändern sich eigentlich die Indices?
Ich bin euch für jede Spur bezügl. der Faltung im Frequenzraum dankbar.
Also das eine Faltung im Ortsbereich einer Multiplikation im Frequenzbereich und umgekehrt entspricht, da sind wir uns einig.
Das zeigt auch folgendes Beispiel:
Code:
x=sin(1:0.1:10); %Signal
h=ones(1,10); % Faltungskern
a=conv(h,x); % Faltung im Ortsraum
b=fft(x,128).*fft(h,128);% Multiplikation im Frequenzraum
c=ifft(b); % Rückführung in den Ortsraum
figure(1) plot(a,'bx-') hold on
plot(real(c),'k')% zu beachten ist, dass wir nur den Realteil betrachten legend('Faltung im Ortsraum','Multiplikation im Frequenzraum')
Ich seh jetzt gerade den Knackpunkt nicht, wo es in deiner Berechnung hakt?
_________________
>> why
The computer did it.
Gast
Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
Verfasst am: 18.04.2009, 12:57
Titel:
Zunächst einmal, möchte ich mich für die rasche Antwort bedanken!
Im Prinzip ähnelt meine Vorgehensweise deinem Beispiel sehr. Allerdings bin ich noch nicht auf die elegante Idee gekommen, die Transformation und die Faltung in einem Schritt durchzuführen. Das Rechteckfenster ist hier zunächst auch nur ein Platzhalter und sollte keine Auswirkungen haben. Später werde ich natürlich ein geeignetes Fenster auswählen müssen. Aber nun gut, das Faltungsergebnis was mir mein Code erzeugt unterscheidet sich von meiner analytischen Lösung. Und die analytische Lösung ist ganz sicher korrekt (da mehrfach abgesichert). Woran könnte das liegen? Ich sehe keine große Abweichung zwischen deinem Beispiel und meinem Code.
Eine Multiplikation im Ortsraum entspricht einer Faltung im Frequenzraum. Du musst also nochmal schauen, wo du jetzt multiplizierst und das bei der analytischen Rechnung bedenken.
Ich hab das mal geplottet und es sieht für mich identisch aus, sowohl mittels conv im ortsraum als auch über die fft.
Zuletzt bearbeitet von Maddy am 18.04.2009, 16:09, insgesamt 2-mal bearbeitet
Gast
Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
Verfasst am: 18.04.2009, 15:11
Titel:
Ich konnte es leider noch nicht prüfen, aber mir ist aufgefallen, dass du beim Plot zweimal das Signal Y1 verwendet hast. Natürlich sehen die Signale dann identisch aus. Ich werde das mal püfen. Vor allem werde vorerst die Fensterung vollständig entfernen.
Arghs, sry das passiert wenn man code verschlimmbessern will.
Ich schau nochmal drüber.
edit: Das "Problem" scheint in der Größe zu liegen. Also dass du quasi ein zu großes Signal über ein anderes Falten möchtest. Die FFT Variante kann aber "quasi" nur den Bereich abdecken, wo die Signale übereinander liegen ohne Überlapp zu haben.
Ich hab das Beispiel oben mal geupdated, wenn du bei h die Größe änderst, wirst du die Abhängigkeit erkennen.
Das die FFT sonstwie krumm aussieht, liegt an numerischen Fehlern (Größenordnung 10^-14 also quasi 0).
_________________
>> why
The computer did it.
Gast
Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
Verfasst am: 18.04.2009, 16:14
Titel:
Hallo,
habe das ganze nun geprüft und ein paar Kleinigkeiten geändert. D. h. die Fensterung habe ich komplett entfernt. Und den überflüssigen Zeitvektor t (tauchte unterhalb deiner Definition nocheinmal auf) entfernt. Weiterhin habe ich Y3 über den Weg deines ersten Beispiels erzeugt. Logischerweise unterscheidet sich das Ergebnis nicht von meiner ursprünglichen Variante, ist aber ohne Frage eleganter. Wenn du dir die Plots anschaust, wirst du feststellen, dass die Faltung im Ortsraum sich von der Faltung im Frequenzbereich unterscheidet. Was könnten wir da noch falsch machen?
Y1FFT = X1FFT .* HFFT;
Y2=conv(x1, h); % normale Faltung im Ortsraum
b=fft(x1,128).*fft(h,128);% Multiplikation im Frequenzraum
% Perform the iFFT
Y1 = real(ifft(Y1FFT)); % <-- hier vorsichtig mit anderen Fenstern, da dies zusätzlich einer Faltung im Ortsraum entspricht.
Y3=real(ifft(b)); % Rückführung in den Ortsraum
Oh, da hat sich mein Beitrag mit deinem überschnitten. Danke für den wichtigen Hinweis mit der Größe werde auch das mal prüfen. (Um Missverständnisse zu vermeiden: Mein Kommentar "...eleganter" bezieht sich natürlich auf deinen Weg, da er ja kürzer ist")
Ist zwar schon eine Weile her als dieses Thema behandelt wurde, aber evtl. interessiert es ja jemanden.
Die Faltung im Frequenzbereich liefert ein anders Ergebnis, da hier eine zyklische Faltung ausgeführt wird. Die Faltung bzw. ihre Inverse sind aber linear. Die Faltung von x[m] und h[n] ergibt einen Ergebnisvektor von y[m+n-1]. Deine beiden Eingangsvektoren sind 100 Werte lang. Die Faltung ergäbe somit 199 Werte. Du führst aber nur eine FFT über 128 Messwerte aus. Die 199 richtigen Werten passen aber nicht in die 128 der FFT .Um also das identische Ergebnis zu Conv() über die Faltung im Frequenzbereich zu erhalten, müssen die Eingangsvektoren x und h schon die Länge m+n-1 haben. In dem letzten Programm genügt es, wenn die FFT über 256 Werte ausgeführt wird.
Es hat hier zwar keinen Einfluss, aber sollte man nicht auch den Imaginärteil für die IFFT benutzen? Das Eingangssignal ist zwar reell, aber eine ungerade Funktion. Somit ist der Imaginärteil nicht null.
---------------------------
Nun habe ich selber aber noch eine Frage. Ich nutze die Faltung u. die inverse Faltung im Frequenzbereich ebenfalls, allerdings in Bezug auf eine Sprungantwort. Mathematisch entsteht eine Sprungantwort aus der Faltung des Sprungs mit der Impulsantwort des Systems. Im Frequenzbereich geht die Faltung dann in die Multiplikation über.
Sprung[f] * Impulsantwort[f] = Sprungantwort[f]
Kenne ich nun den Eingang auf das System und die Systemantwort, kann man durch Divison die Frequenzantwort berechnen.
Frequenzantwort[f] = Sprungantwort[f] / Sprung[f]
Durch IFFT erhalte ich dann meine Impulsantwort. Mit FFT Faltung bekomme ich hier leider noch kein Ergebnis, da Division durch 0. Die Funktion deconv() liefert mir aberGenauso kann ich auch den Eingang bereichnen, wenn ich die anderen beiden Teile habe.
Sprung[f] = Sprungantwort[f] / Frequenzantwort[f]
Dies gelingt mit der FFT Faltung, aber nicht mit [q,r]=deconv(). Der Quotient ergibt dann 1 und der Rest ist das gezeigte Signal. Außerdem ergeben sich Unterschiedliche Ergebnisse für den s- und z-Bereich.
Code:
N = 1024;
Ta = 1/613;
Auswahl = 2;
%Zeitvektor
T = 0:Ta:((N-1)*Ta) ;
%Streckenparameter
D1=0.1;
T1=0.01;
Kst=1;
%zunächst im Bildbereich
s = tf('s');
Fst = (Kst/(T1^2*s^2 + 2*D1*T1*s + 1));
%z-Transformation
z = tf('z',Ta);
Fzt = c2d(Fst,Ta);
for i=1:(length(x1)),
if(x1(i)== 0) %Division durch 0 vermeiden
x1(i)= 0.000001;
end end for i=1:(length(h1)),
if(h1(i)== 0) %Division durch 0 vermeiden
h1(i)= 0.000001;
end end
% Perform the FFT
X1FFT = fft(x1,2048);
HFFT = fft(h1,2048);
%Auswahl: 1 = Faltung, 2 = inverse Faltung
Auswahl=2;
switch(Auswahl) case1 % Perform the convolution
Y1FFT = X1FFT .* HFFT;
Y2=conv(x1, h1); % normale Faltung im Ortsraum
b=fft(x1,2048).*fft(h1,2048);% Multiplikation im Frequenzraum case2 % Perform the convolution
Y1FFT = X1FFT ./ HFFT;
[q,Y2]=deconv(x1, h1); % normale Faltung im Ortsraum
b=fft(x1,2048)./fft(h1,2048);% Multiplikation im Frequenzraum end
% Perform the iFFT
Y1 =real(ifft(Y1FFT)); % <-- hier vorsichtig mit anderen Fenstern, da dies zusätzlich einer Faltung im Ortsraum entspricht.
Y3=(ifft(b,'symmetric')); % Rückführung in den Ortsraum
Du kannst Beiträge in dieses Forum schreiben. Du kannst auf Beiträge in diesem Forum antworten. Du kannst deine Beiträge in diesem Forum nicht bearbeiten. Du kannst deine Beiträge in diesem Forum nicht löschen. Du kannst an Umfragen in diesem Forum nicht mitmachen. Du kannst Dateien in diesem Forum posten Du kannst Dateien in diesem Forum herunterladen
MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.
.Um also das identische Ergebnis zu Conv() über die Faltung im Frequenzbereich zu erhalten, müssen die Eingangsvektoren x und h schon die Länge m+n-1 haben. In dem letzten Programm genügt es, wenn die FFT über 256 Werte ausgeführt wird.
RSS