Hilberttransformation: Filtern&Rückgewinnung des Signal
Michi8888
Gast
Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
Verfasst am: 27.09.2008, 19:00
Titel: Hilberttransformation: Filtern&Rückgewinnung des Signal
Ich möchte eine Funktion x mit dem Hilbertfilter fh filtern und anschließend daraus wieder das Signal x rückgewinnen?
Das Filtern hätte ich so gemacht:
"y= filter(fh,1,d);"
Stimmt das?
Und wie kann ich nun aus y wieder das ursprüngliche Singal x gewinnen?
um die Faltung bzw die Filterung rückgängig zu machen, kannst du y und fh per fft in den Frequenzbereich bringen und dann dividieren. Das entspricht dem inversen der Faltung (Deconvolution).
Matlab besitzt auch noch eine Funktion "deconv", diese basiert glaube ich auf einer Polynomdivision, ich habe aber noch nicht damit gearbeitet. Kannst ja mal schauen. Stichwort ist auf jeden Fall "deconvolution".
mfg
Carsten
Michi8888
Gast
Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
Verfasst am: 30.09.2008, 10:19
Titel:
Vielen vielen Dank Carsten!
So wie du es geschrieben hast ist es einleuchtend für kontinuierliche Signale. Leider brauche ich es für zeitdiskrete Signale.
Hier mein Quellcode (wichtig ist eigentlich nur der Teil in der Mitte! Damit man es in Matlab ablaufen lassen kann hab ich alles angegeben.):
"M = 100;
n = -M:M; % Zeitachse
tg = 1/2; % Grenzfrequenz angegeben in 1/pi
Nf = 512; % Anzahl der Frequenzkomponenten, da ich später die fft auf zeitdiskrete Signale anwende
fr = linspace(-1,1,Nf); % Frequenzachse in theta/pi
f = tg*sinc(tg*n); % nur eine Funktion die ich später verwende
x1 = tg*sinc(tg*n/2); %auch nur eine Funktion die ich später verwende
fh = f .*x; % mit dieser Funktion möchte ich filtern
d = (x1 .* x1) .* exp(j*pi*n/2); % dieses Signal soll gefiltert werden
dd = filter(fh,1,d); % Filterung
D = fft(d,Nf);
DD = fft(dd,Nf);
FH = fft(fh,Nf);
DDR = FH./DD;
ddr = ifft(DD, Nf); % Denke, dass hier der Fehler liegt.
%%% alles ab hier dient nur der Darstellung
subplot(3,2,1), plot(n, d);
title('Signal d');
xlabel('n'); ylabel('d');
x1 gemeint war, ansonsten kann ich deinen code nicht ausführen . Der Fehler wird bei mir von der plot-Funktion zurückgegeben.
Zitat:
??? Error using ==> plot
Vectors must be the same lengths.
Error in ==> forumtest at 39
subplot(3,2,5), plot(n, ddr);
Das bedeutet, dass n und ddr nicht gleich lang sind. Das müssen sie aber sein da er ja ddr über n plotten soll. Du musst den vektor ddr beschneiden, damit das passt. Du hast ihn ja vorher auf Nf=512 samples aufgeblasen. Die fft macht das mit zeropadding, also kannst du theretisch einfach sagen
Danke! Stimmt, damit is die Fehlermeldung weg. Ich habe die Zeile
Zitat:
x = x1 .*x1;
vergessen, daher ließ es sich nicht ausführen. (In einer Plotfunktion hatte ich auch einen kleinen Fehler; den habe ich jetzt beseitigt.) Unabhängig davon sollte dann ddr wieder gleich dem Signal vor der Filterung, also gleich d sein. Leider ist dem nicht so und ich verstehe nicht warum bzw. was noch fehlt.
Code:
M = 100;
n = -M:M; % Zeitachse
tg = 1/2; % Grenzfrequenz angegeben in 1/pi
Nf = 512; % Anzahl der Frequenzkomponenten
fr = linspace(-1,1,Nf); % Frequenzachse in theta/pi
f = tg*sinc(tg*n); % Zeitsignal f
x1 = tg*sinc(tg*n/2);
x = x1 .*x1;
fh = f .*x; % Multiplikation zwecks Phasenverschiebung
d = x .*exp(j*pi*n/2);
D = fft(d,Nf);
dd = filter(fh,1,d);
DD = fft(dd,Nf);
FH = fft(fh,Nf);
DDR = FH./DD;
ddr = ifft(DDR, Nf); % ??
Hi,
habs mir mal angeschaut. Also zum einen stimmt was nicht der Funktion filter, die du verwendest. Diese liefert dir nur das halbe Ergebnis der Faltung. Bei der Faltung zweier Vektoren a 201 samples sollte man einen Vektor bekommen der 401 samples lang ist (Einlaufzeit+1 Sample wo sie übereinanderliegen+Auslaufzeit). Du kannst die Vektoren mit der Funktion "conv" direkt falten und dir dann die mittleren 201 samples ausschneiden, um das Ergebnis darzustellen. Oder du verlängerst deinen Zeitvektor, denn eine Filterung hat nunmal eine Verzögerung aufgrund der Schieberegister (mal in Digitaltechnik
gedacht ).
Der Grund warum du bei der Deconvolution nicht das richtige heraus bekommst ist, dass du FH./DD rechnest. Das macht ja keinen Sinn, denn du willst ja nicht dein Faltungsergebnis aus dem Filter "rausfalten", sondern umgekehrt. Dann noch zur Darstellung der Spektren, dabei musst du beachten, dass die Funktion fft dir die Frequenzen von 0.1.2....Nyquist...-2.-1 liefert (Das hat damit zu tun, dass die DFT periodische Spektren liefert). Die Funktion fftshift klappt den negativen Frequenzanteil dahin wo er hingehört. Ich hab deinen code ein wenig verändert und dir noch ein zusätzliches plot-fenster gemacht.
Code:
M = 100;
n = -M:M; % Zeitachse
tg = 1/2; % Grenzfrequenz angegeben in 1/pi
Nf = 512; % Anzahl der Frequenzkomponenten
fr = linspace(-1,1,Nf); % Frequenzachse in theta/pi
f = tg*sinc(tg*n); % Zeitsignal f
x1 = tg*sinc(tg*n/2);
x = x1 .*x1;
fh = f .*x; % Multiplikation zwecks Phasenverschiebung
d = x .*exp(j*pi*n/2);
D = fft(d,Nf);
%dd = filter(fh,1,d); Diese Funktion liefert nur das halbe Ergebnis der
%Faltung
dd=conv(fh,d); %Dieser Vektor ist 401 samples lang (Faltung=>length(d)+length(fh)-1)
DD = fft(dd,Nf);%Die fft füllt dd (401 samples) auf 512 samples auf (zeropadding)
FH = fft(fh,Nf);%fh (201 samples) -> 512 samples
DDR = DD./FH; %nicht FH./DD
ddr = ifft(DDR, Nf); % 512 Spektralanteile => 512 samples im Zeitbereich
Danke!!!
Eine Frage (die sich mir auch bei anderen Aufgaben gestellt hat) habe ich noch: Warum darf man fftshift anwenden? fftshift entspricht doch einer Phasenverschiebung und somit einer Modifikation des Spektrums oder?
lg Michi
also von einer Phasenverschiebung weiss ich nichts. Die matlab-Hilfe sagt folgendes:
Zitat:
FFTSHIFT Shift zero-frequency component to center of spectrum.
For vectors, FFTSHIFT(X) swaps the left and right halves of
X. For matrices, FFTSHIFT(X) swaps the first and third
quadrants and the second and fourth quadrants. For N-D
arrays, FFTSHIFT(X) swaps "half-spaces" of X along each
dimension.
FFTSHIFT(X,DIM) applies the FFTSHIFT operation along the
dimension DIM.
FFTSHIFT is useful for visualizing the Fourier transform with
the zero-frequency component in the middle of the spectrum.
Nach meinem Verständnis wird die Phasenlage einzelner Frequenzanteile dadurch nicht beeinflusst. Man ändert ja nur Reihenfolge der Anteile nicht aber dessen Phase. Worauf man achten muss ist, dass man das geshiftete Spektrum nicht ohne weiteres per ifft zurücktransformieren darf. Dann kommt es zu einer Verfälschung des Signals im Zeitbereich. Für die Darstellung in einer Plot-Funktion ist es jedoch anschaulicher.
mfg
Carsten
Michi8888
Gast
Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
Verfasst am: 01.10.2008, 14:07
Titel:
Naja, die Matlab-Funktionsbeschreibung habe ich eh auch schon gelesen aber ist mir nicht klar wie ich das in folgendem Bsp. deuten soll.
Beim folgenden Bsp. soll ich veranschaulichen was bei Multiplikation der Zeitsignale f[n] und x[n] geschieht. Diese Aufgabe wird im Zuge einer Ausarbeitung über die Hilberttransformation behandelt, daher nehme ich an, dass das Ergebnis das Hilbert-Filter, welches 1 für 0<theta<pi und 0 für -pi<theta<0 ist.
Das Ergebnis ist das Signal mod, dessen Spektrum MOD laut meiner Berechnung dem eines Hilbertfilters gleicht - nur mit dem Unterschied, dass es negative Frequenzkomponenten aufweißt. Wende ich auf MOD die Funktion fftshift an, so wandert das Spektrum nach rechts und das Ergebnis ist das Hilbertfilter. ... Kurz gesagt: Es hängt von fftshift ab, ob das Ergebnis das Hilbertfilter ist oder nicht.
Hier der Quellcode, wobei das letzte geplottete Spektrum das Ergebnis darstellt. (Alle anderen Funktionen sind denke ich irrelevant.)
Zitat:
M = 100;
n = -M:M; % Zeitachse
tg = 1/2; % Grenzfrequenz angegeben in 1/pi
f = tg*sinc(tg*n); % Zeitsignal f
subplot(3,2,1), plot(n, f);
title('Zeitsignal f');
xlabel('n'), ylabel('f'), grid on;
Nf = 512; % Anzahl der diskreten Fourierkoeffizienten
F = fftshift(fft(f,Nf)); % DFT und Zentrierung des Spektrums
fr = linspace(-1,1,Nf); % Frequenzachse in theta/pi
subplot(3,2,2), plot(fr,abs(F(1:Nf)));
title('Spektrum F von f');
xlabel('\theta/\pi'), ylabel('|F(e^{j\theta})|'), grid on;
x = exp(j*pi*n/2); % Phasenverschiebungsterm
n1 = 0:.1:3*pi; % n1, x1 werden nur verwendet, um den Graphen
% der Funktion x anschaulich zu gestalten
x1 = exp(j*pi*n1/2);
subplot(3,2,3), plot3(real(x1), imag(x1), n1);
title('Signal x für die Phasenverschiebung');
xlabel('real(x)'), ylabel('imag(x)'), zlabel('n'), grid on;
% subplot(3,2,3), plot3(real(x), imag(x), n), xlabel('real(x)'),
% ylabel('imag(x)'), zlabel('n'), grid on;
% Wäre die korrekte Version,ist aber nicht anschaulich, da n sehr groß ist.
X = fftshift(fft(x,Nf)); % DFT und Zentrierung des Spektrums
subplot(3,2,4), plot(fr,abs(X(1:Nf)));
title('Spektrum X von x');
xlabel('\theta/\pi'), ylabel('|X(e^{j\theta})|'), grid on;
mod = f .*x; % Modulation des Signals f mit dem Phasenver-
% verschiebungsterm
subplot(3,2,5), plot(n, mod), xlabel('n'),
title('Moduliertes Zeitsignal');
ylabel('mod'), grid on;
MOD = fftshift(fft(mod, Nf)); %%%%%% DFT und Zentrierung des Spektrums - Wenn ich fftshift nicht anwende hat das Spektrum
negative statt positive Spektralanteile
subplot(3,2,6), plot(fr,abs(MOD(1:Nf)));
xlabel('\theta/\pi'), ylabel('|MOD(e^{j\theta})|'), grid on;
Falls du Zeit hast dir das noch anzusehen, würd ich mich natürlich sehr freuen. ... Denn morgen habe ich Abgabe.
lg Michi
Michi8888
Gast
Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
Verfasst am: 01.10.2008, 14:11
Titel:
Hier nochmals derselbe Code. Habe ihn vorhin versehentlich als Zitat eingefügt.
Code:
M = 100;
n = -M:M; % Zeitachse
tg = 1/2; % Grenzfrequenz angegeben in 1/pi
f = tg*sinc(tg*n); % Zeitsignal f subplot(3,2,1), plot(n, f);
title('Zeitsignal f');
xlabel('n'), ylabel('f'), grid on;
Nf = 512; % Anzahl der diskreten Fourierkoeffizienten
F = fftshift(fft(f,Nf)); % DFT und Zentrierung des Spektrums
fr = linspace(-1,1,Nf); % Frequenzachse in theta/pi subplot(3,2,2), plot(fr,abs(F(1:Nf)));
title('Spektrum F von f');
xlabel('\theta/\pi'), ylabel('|F(e^{j\theta})|'), grid on;
x = exp(j*pi*n/2); % Phasenverschiebungsterm
n1 = 0:.1:3*pi; % n1, x1 werden nur verwendet, um den Graphen % der Funktion x anschaulich zu gestalten
x1 = exp(j*pi*n1/2);
subplot(3,2,3), plot3(real(x1), imag(x1), n1);
title('Signal x für die Phasenverschiebung');
xlabel('real(x)'), ylabel('imag(x)'), zlabel('n'), grid on;
% subplot(3,2,3), plot3(real(x), imag(x), n), xlabel('real(x)'), % ylabel('imag(x)'), zlabel('n'), grid on; % Wäre die korrekte Version,ist aber nicht anschaulich, da n sehr groß ist.
X = fftshift(fft(x,Nf)); % DFT und Zentrierung des Spektrums subplot(3,2,4), plot(fr,abs(X(1:Nf)));
title('Spektrum X von x');
xlabel('\theta/\pi'), ylabel('|X(e^{j\theta})|'), grid on;
mod = f .*x; % Modulation des Signals f mit dem Phasenver- % verschiebungsterm subplot(3,2,5), plot(n, mod), xlabel('n'),
title('Moduliertes Zeitsignal');
ylabel('mod'), grid on;
MOD = fftshift(fft(mod, Nf)); % DFT und Zentrierung des Spektrums, wenn ich fftshift nicht anwende habe ich negative statt positive Spektralanteile subplot(3,2,6), plot(fr,abs(MOD(1:Nf)));
xlabel('\theta/\pi'), ylabel('|MOD(e^{j\theta})|'), grid on;
ehrlich gesagt verstehe ich dein Problem nicht so ganz. Es kommt doch das raus, was du erwartest, wenn du fftshift anwendest. Am Spektrum von f siehst du auch, das fftshift notwendig ist, um das Spektrum korrekt darzustellen, denn die Fouriertransformierte der sinc-Funktion ist ein Rechteck und umgekehrt. Würdest du bei der Darstellung von F nicht fftshift anwenden, würdest du zwei Hälften eines Rechtecks sehen, was nicht richtig ist (zumindest bei deiner Achsen-Skalierung). Die fft erzeugt kein Spektrum, so wie man es in Büchern sieht. Das Spektrum der DFT ist periodisch (in der Theorie). Die fft liefert dir praktisch zuerst die zweite Hälfte der ersten Periode und dann die erste Hälfte der zweiten Periode (Ich denke das liegt daran, dass Matlab keine negativen Indizes kennt). Deshalb ist die Funktion fftshift notwendig, um ein Spektrum zu erhalten, was dem aus Büchern/Vorlesung entspricht.
mfg
Carsten
Michi8888
Gast
Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
Verfasst am: 01.10.2008, 19:15
Titel:
Genau das wars, was mir nicht klar war. Jetzt habe ich es gecheckt! Ich habe nicht daran gedacht, dass es einen Unterschied zwischen fft und der DFT gibt.
Vielen Dank nochmals für deine Hilfe!!
lg Michi
Einstellungen und Berechtigungen
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.
. Der Fehler wird bei mir von der plot-Funktion zurückgegeben.
RSS