Da ich nicht gerade ein Profi bin, wenn es um Signalverarbeitung geht, muss ich mich mal wieder an euch wenden.
Hintergrund:
Ich habe ein Programm geschrieben, mit dem ich gemessene Zeitdaten (Audiodaten) auswerten kann. Ich lese also die Zeitdaten ein und führe das Programm aus. Das funktioniert auch wunderbar. Die Ergebnisse sind richtig.
Nun möchte ich das Programm gerne auf selbst generierte Daten anwenden. Dazu liegt mir die entsprechende Übertragungsfunktion (G) vor (im Frequenzbereich).
Ich gehe also wie folgt vor:
1. Generieren von Signal x (Zeitbereich)
2. mittels fft transformieren --> Signal X (Frequenzbereich)
3. Multiplikation von G und X --> Ausgangssignal Y (Frequenzbereich)
4. mittels ifft transformieren --> Ausgangssignal y (Zeitbereich)
5. mit y das Programm ausführen
Nur leider passt das, was bei dem Programm dann rauskommt überhaupt nicht. Und da das Programm ja mit den "richtigen" Daten einwandfrei funktioniert, gehe ich davon aus, dass ich bei der Generierung von y etwas falsch mache.
% Signal erzeugen (Zeitbereich)
x=-1+2.*rand(1,length(t));
hf=design(fdesign.bandpass... % Filterparameter für ('N,F3dB1,F3dB2',100,...% ein Bandpassfilter zwischen 20,12000,fs)); % 20 Hz und 12 kHz
x=filter(hf,x);
% Darstellung des Zeitsignals figure(1) subplot(4,1,1) plot(t,x) title('Eingangssignal (Zeitbereich)')
% mittels FFT in den Frequenzbereich transformieren
X=fft(x,nfft);
X_abs=abs(X);
% Darstellung des Spektrums subplot(4,1,2) plot(fs/nfft*(0:(nfft/2-1)),X_abs(1:nfft/2)) title('Spektrum des Eingangssignals')
% Multiplikation des Spektrums mit der Übertragungsfunktion
Y=G*X;
Y_abs=abs(Y);
% Darstellung des veränderten Spektrums subplot(4,1,3) plot(fs/nfft*(0:(nfft/2-1)),Y_abs(1:nfft/2)) title('Spektrum des Ausgangssignals')
% mittels IFFT in den Zeitbereich zurücktransformieren
y=ifft(Y,nfft);
% Darstellung des wiederhergestellten Zeitsignals subplot(4,1,4) plot(t,y(1:length(t))) title('Ausgangssignal (Zeitbereich)')
Nun zu den konkreten Fragen:
I.) Sind die Schritte 1.-5. theoretisch richtig?
II.) Muss das Rauschen (x) gefiltert werden oder gibt es eine andere Lösung ein bandbegrenztes Rauschen zu erzeugen?
III.) Ist das Spektrum des Einganssignals (X) richtig oder müsste ich es wegen der Spiegelung eigentlich noch mit 2 multiplizieren?
IV.) Ist es korrekt die Daten nach der ifft einfach abzuschneiden und den Rest nicht zu beachten? D.h. die hier geplotteten Daten möchte ich als Zeitdaten mit dem Programm auswerten, ist das richtig?
Ich hoffe das war einigermaßen verständlich. Ansonsten bitte einfach nochmal nachfragen.
nach der ifft habe ich doch gar kein abs() gemacht.
Der real() Befehl nach der ifft hat aber nur eine Auswirkung, wenn ich eine komplexe Übertragungsfunktion habe, oder?
Auf jeden Fall ein guter Tipp. Weil die habe ich eigentlich auch.
Klärt nur leider nicht die eigentlichen Fragen und löst das Problem nicht. Ich bekomme damit noch immer keine sinnvollen Ergebnisse.
ich habe dir mal dein Programm etwas umgeändert und gleich noch die richtige Skalierung der Amplitude im Frequenzbereich eingebaut. Das Programm liefert in y ein vollkommen plausibles Ergebnis. Ich habe der Einfachheit halber einen sinus für das Signal x gewählt. Der Bandpass macht hier so natürlich wenig Sinn...aber darauf kommt es erstmal nicht. Die Multiplikation im Frequenzbereich entspricht ja einer Faltung im Zeitbereich. Mit G = 0.5 wird die Amplitude von x halbiert. Genau dieses Ergebnis ist nach der ifft zu sehen.
% Signal erzeugen (Zeitbereich) % x=-1+2.*rand(1,length(t));
x = sin(2*pi*100*t); % für einfachen Test ein Sinus mit 100 Hz
hf=fdesign.bandpass... % Filterparameter für ('N,fc1,fc2',100,...% ein Bandpassfilter zwischen 20,12000,fs); % 20 Hz und 12 kHz
d = butter(hf);
x=filter(d,x);
% Darstellung des Zeitsignals figure(1) subplot(2,1,1) plot(t,x)
%title('Eingangssignal (Zeitbereich)')
% mittels FFT in den Frequenzbereich transformieren
X=fft(x,nfft);
k = nfft/2 + 1; % Anzahl diskrete Frequenzstellen im pos. Bereich
df = fs/nfft; % Spektrumauflösung
f = 0:df:fs/2; % Frequenzvektor für positiven Freq.bereich
% Skalierung des Betrag % Gleichsignalanteil und Amplitude bei fs/2 müssen extra skaliert werden
magY = [Fy_pos0(1)/nfft ;Fy_pos0(2:k-1)/(nfft/2);Fy_pos0(k)/nfft];
% Umrechnung in dB
magY_dB = 20*log10(magY + eps); % eps = kleine Konstante zur Vermeidung von log(0)
% Darstellung des veränderten Spektrums subplot(2,1,2) hold on;
plot(f,magY_dB,'r') grid on;
%title('Spektrum des Ausgangssignals') legend('Eingang','Ausgang');
% mittels IFFT in den Zeitbereich zurücktransformieren
y=ifft(Y,nfft);
% Darstellung des wiederhergestellten Zeitsignals subplot(2,1,1) hold on;
plot(t,y(1:length(t)),'r.-') grid on;
title('Ein- und Ausgangssignal (Zeitbereich)') legend('Eingang','Ausgang')
Die Faltung funktioniert so aber nur in diesem trivialen Bsp. für G. Bei einer komplexen Funktion wohl nicht mehr. Normalerweise würde man die Impulsantwort des Systems nehmen und diese neben dem Signal auch in den Frequenzbereich transformieren. Dabei muss vorher die Ausgangslänge eines Signals bei der Faltung berücksichtigt werden.
x[n] = Eingangssignal mit n Werten
h[m] = Impulsantwort des Systems mit m Werten
y[n + m - 1] = x[n] * h[m] = Länge des Ausgangssignal.
x und h MÜSSEN unbedingt noch vor der FFT auf die Länge n+m-1 durch Anhängen von Nullen gebracht werden. Hier siehst du so eine Faltungsfunktion im Freq. bereich.
Code:
function output = FFT_Faltung(sig1, sig2) % Die Funktion führt eine schnelle Faltung mittels FFT aus % Input: % sig1 = Messsignal % sig2 = Filterimpulsantwort % Output: % output = gefiltertes Messignal
%--------------------------------------------------------------------------
% kopieren der Eingangsvektoren
sig1 = double(sig1(:));
sig2 = double(sig2(:));
% Faltungssatz sig(m) * h(n) = output(m+n-1):
outlength = length(sig1)+length(sig2)-1;
% nächste Zweierpotenz damit FFT Algorithmus verwendet wird
fftsize = 2^nextpow2(outlength);
% Messsignal % mit Nullen auf fftsize auffüllen
sig1 = [sig1; zeros(fftsize-length(sig1),1)];
% Diese Zeile könnte man sich auch sparen, da die fft(sig1,fftsize) schon automatisch mit Nullen auffüllt, wenn sig1 < fftsize ist. Die Zeile ist nur zur Verdeutlichung was Zeropadding ist.
%Berechnung des Frequenzspektrums
sig1 = fft(sig1,fftsize);
% Filter - Impulsantwort % mit Nullen auf fftsize auffüllen
sig2 = [sig2; zeros(fftsize-length(sig2),1)];
% Berechnung des Frequenzspektrums
sig2 = fft(sig2,fftsize);
% Faltung durch Multiplikation der Spektren % Rücktransformation in den Zeitbereich
conv_raw = ifft(sig1.*sig2);
% Ausgabe
output = transpose(conv_raw(1:outlength));
Hallo,
jetzt habe ich doch nochmal eine Frage zu dem Thema.
Und zwar geht es genau um das, was Andreas schonmal angesprochen hat. Nämlich den Umgang mit den Daten nach der IFFT.
Wenn ich an dieser Stelle den abs()-Befehl verwende, dann erhalte ich die erwarteten Ergebnisse. Mit dem real()-Befehl stimmen die Ergebnisse nicht.
Laut der sogenannten "Zeitkonvention" müsste ich aber den real()-Befehl verwenden um aus den komplexen Amplituden wieder Orts- und Zeitverläufe zu erhalten:
Kann mir jemand erklären, warum der abs()-Befehl hier die richtigen Ergebnisse liefert?
Hi,
in dem geposteten Code verwende ich weder abs() noch real(). Richtig.
Aber wie beschrieben, ist die eigentliche Übertragungsfunktion komplizierter. Sie ist beinhaltet die e-Funktion und ist komplex. Und mein y ist dann jedenfalls auch komplex.
Sorry, wenn das missverständlich war.
Ich versuche die Ergebnisse aus einem Paper nach zu vollziehen und programmiere das in Matlab nach (naja - ich versuche es zumindest). Und wenn ich nach der IFFT den abs()-Befehl verwende, dann sehen meine Ergebnisse so aus wie in dem Paper.
Ich vermute, dass es sich hier bei mir einfach um ein Signalverarbeitungs-Verständnisproblem handelt, oder?
Die Multiplikation im Frequenzraum entspricht ja einer Faltung. Führt denn die Faltung von G und x im Zeitbereich auch zu einem komplexen Signal? Wenn ja, darf sich doch auch durch den Umweg über den Frequenzbereich daran nichts ändern.
Poste mir doch mal bitte ein Bsp. für deine Übertragungsfunktion.
Hi,
jetzt weiß ich nicht genau, ob ich dich richtig verstanden habe...
Also wenn ich das (Zeit-) Signal x mit der Impulsantwort g (IFFT(G)) falte, dann ist das Ergebnis auch komplex. Genau, wie bei der Multiplikation im Frequenzbereich auch.
Aber wie hilft mir diese Erkenntnis nun weiter?
Die Übertragungsfunktion sieht z.B. so aus:
Dabei sind die x Skalare, k die Wellenzahl und j die imaginäre Einheit.
% Signal erzeugen (Zeitbereich)
x = sin(2*pi*10*t);
% mittels FFT in den Frequenzbereich transformieren
nfft=2^nextpow2(length(x));
X=fft(x,nfft);
% Frequenzvektor für die frequenzabhängige Übertragungsfunktion bestimmen
df=fs/nfft; % Auflösung des Spektrums
f=0:df:fs-df; % Vektor mit length(f)=nfft
Das Ausgangssignal y im Zeitbereich ist komplex.
Meine Messdaten sind aber reell und ich möchte auch reelle generierte Daten haben.
Was ist nun der richtige Weg um aus den komplexen Daten reelle Daten zu machen?
ifft tests X to see whether vectors in X along the active dimension are conjugate symmetric. If so, the computation is faster and the output is real. An N-element vector x is conjugate symmetric if x(i) = conj(x(mod(N-i+1,N)+1)) for each element of x.
Das Signal X(f) hat nach der FFT diese Form. Hier ist nämlich von 0...nfft/2 + 1 der positive Frequenzbereich (0...+fs/2 [Hz]) und von nfft/2+1...N der negative Frequenzbereich (-fs/2...0-df [Hz]) enthalten, wobei die Frequenz 0 und fs/2 in dem Vektor X nur einmal vorhanden sind. Der negative Bereich ist aber einfach eine Spiegelung (konjugiert komplex) und die Amplitude somit nur die Hälfte, weshalb man anders Skalieren muss, will man nur den pos. Freq.-raum in einem Graph darstellen
Du hast nun aber ein Frequenzantwort von G(f) im Bereich 0...Fs erstellt und das kann nicht passen.
Code:
% Frequenzvektor für die frequenzabhängige Übertragungsfunktion bestimmen
df=fs/nfft; % Auflösung des Spektrums % !!! Änderung auf fs/2 !!!
f1=0:df:fs/2; % Vektor mit length(f)=nfft
% Multiplikation des Spektrums mit der Übertragungsfunktion
l = nfft/2 + 1;
% es wird von X nur der Bereich 0...fs/2 [Hz] verwendet
Y1=transpose(G.*X(1:l));
% spektral spiegeln
Y=[Y1(1:l); conj(flipud(Y1(2:l-1)))];
Allerdings erhalte ich so immernoch ein komplexes Signal y. Ich muss da auch mal noch weiter grübeln. Ist es richtig, dass die Funktion G einen Gleichrichter darstellt?
% mittels FFT in den Frequenzbereich transformieren
nfft=2^nextpow2(length(x));
X=fft(x,nfft);
% Frequenzvektor für die frequenzabhängige Übertragungsfunktion bestimmen
df=fs/nfft; % Auflösung des Spektrums
f=-fs/2:df:fs/2-df; % Frequenzvektor mit length(f)=nfft
f=ifftshift(f); % Umordnung der Frequenzen, damit die % Reihenfolge zu den Werten der fft % passt
Wenn jetzt noch jemand erklären kann, warum y komplex ist, wäre das super. Bei einem Rauschen als Eingang wird der komplexe Anteil übrigens noch größer.
G ist übrigens kein Gleichrichter. Es geht um ein Mikro.
Zuletzt bearbeitet von laupl am 25.04.2013, 15:49, insgesamt einmal bearbeitet
Wenn man abs(y) darstellt, werden die negativen Halbwellen des Sinus umgekehrt. Du sprachst ja mal davon, dass dein Ergebnis mit abs(y) deinen Erwartungen entspricht. Was ich dann so sehe, ist ein Zweiweg-Gleichrichter . Mit real(y) ergibt sich natürlich etwas anderes.
Das Problem liegt nach wie vor daran, dass Y nicht den geforderten Aufbau hat, wie in dem Befehl zur ifft beschrieben wird.
Zitat:
An N-element vector x is conjugate symmetric if x(i) = conj(x(mod(N-i+1,N)+1)) for each element of x.
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.
. Mit real(y) ergibt sich natürlich etwas anderes.
RSS