Verfasst am: 22.05.2012, 15:12
Titel: Filterung von Signalen durch tfestimate
Hi,
nachdem ich nun einige Zeit im Internet nach einer geeigneten Lösung gesucht habe, aber leider nicht fündig geworden bin, wollte ich hier mal nachfragen.
Ich habe ein Eingangssignal und ein Ausgangssignal, welches durch die umgebende Struktur unterschiedlich gedämpft bzw. verstärkt wird und möchte die Übertragungsfunktion der Struktur bestimmen, um diese später auf weitere Eingangssignale anwenden zu können, um somit deren Ausgangssignal, ohne diese in die umgebende Struktur einbauen zu müssen, "voraussagen" zu können.
Nach einiger Suche bin ich auf den Befehl tfestimate gestoßen, mit dem ich die Übertragungsfunktion schätzen kann.
Nun stehe ich aber vor dem Problem, wie ich den durch tfestimate erhaltenen Vektor einsetzen muss, um das Signal zu filtern.
Von dem filter-Befehl bin ich es gewohnt zwei Vektoren a und b anzugeben. Durch tfestimate erhalte ich nur einen Vektor.
Irgendwo hatte ich gelesen dass eine Filterung über
aber wenn man sich die FFT des so erzeugten Ausgangssignals anschaut, stimmen die Verstärkungen und Dämpfungen nicht mit denen der Übertragungsfunktion überein.
Kann mir jemand von euch sagen, wie ich also die Werte von tfestimate verwenden muss, um ein Signal damit filtern zu können?
Vielen Dank für eure Hilfe schon mal im Voraus!
Thomas
tfestimate
Gast
Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
Verfasst am: 23.05.2012, 10:41
Titel:
Ich bins nochmal:
Gibt es denn eine Möglichkeit oder einen Befehl, aus zwei Signalen, direkt die Filterparameter bestimmen zu können?
Die Übertragungsfunktion benötige ich nicht unbedingt.
Ziel ist es, wie bereits oben erwähnt, aus einem gegebenen Eingangs- und Ausgangssignal einen Filter zu erzeugen, mit dem auch weitere Eingangssignale gefiltert werden können, um so die Ausgangssignale zu erhalten.
Hier gibt es einige Möglichkeiten...ich gehe mal auf 2 ein:
1.) Das Ausgangssignal y[n] ist ja das Ergebnis der Faltung von Eingangssignal x[n] mit der Impulsantwort h[m] deines Systems. Man kann diese "Verkettung" nutzen und sozusagen eine "Enthaltung" durchführen. Dabei nutzt man einen Vorteil der durch die Transformation in den Freq.raum entsteht. Hier wird die Faltung zu einer einfachen Multiplikation.
x[n] * h[m] = y[n+m-1]
h[m] = y[n+m-1] /x[n]
Du transformierst also y und x mittels FFT getrennt in den Freq.-bereich. Achte dabei auf die Länge n+m-1, es müssen beide Signale gleich lang sein. Falls nicht vorher Nullen anhängen. Dann führst du die Division durch und erhältst die Frequenzantwort H[f] deines Filters. Mit der IFFT bekommst du dann die Impulsantwort. Oder aber du transformierst zukünftig alle Eingangssignale in den Freq.raum und machst die Faltung dort. Das ist ohnehin schneller...das Ergebnis (=Ausgangssignal) mittels IFFT wieder in den Zeitbereich.
2.) Mittels eines ARX und ARMAX die Parameter des Systems schätzen. Es gibt in Matlab auch noch weitere Modelle. Ich habe zu dem Thema auch schon einiges gepostet, wie man da vorgeht.
z.B. hier: http://www.gomatlab.de/frage-zu-arx.....9877,highlight,armax.html
Allerdings muss man vorher die Ordnung des Systems festlegen und hier braucht es entweder Erfahrung oder viel Zeit um die beste Approxm. zu finden.
tfestimate
Gast
Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
Verfasst am: 23.05.2012, 16:32
Titel:
Hi DSP,
vielen Dank für deine Antwort und die guten Links.
Variante 2 scheidet bei mir aus, da ich die entsprechende Toolbox nicht habe, aber mit Variante 1 funktioniert es perfekt. Das letzte mal als ich Variante 1 getestet hatte war es irgendwie an der Rücktransformation gescheitert, da ich dort keine sinnvollen Ergebnisse erhalten hatte, aber nun klappt es mit dem Code
Code:
%Berechnung der Frequenzantwort
H=fft(out)./fft(in); %out und in sind beide gleich lang
%zu Testzwecken erneute Berechnung von out durch Multiplikation von Eingangssignal und Frequenzantwort im Frequenzbereich
OUTTEST=fft(in).*H;
%Umrechnung in Zeitbereich
outtest=ifft(OUTTEST);
Für meinen Zweck würde ich die Frequenzantwort gerne hinterlegen, um später weitere Ausgangssignale berechnen zu können. Ich habe das Ganze nun mal für ein längeres Eingangssignal getestet und H einfach mit Nullen auf die entsprechende Länge gebracht. Das Ergebnis weicht allerdings deutlich von den Erwartungen ab.
Wie genau gehe ich hierbei vor?
Was ich prinzipiell ja möchte, ist dass ich zu jeder Frequenz (bis 24kHz) eine Verstärkung hinterlegt habe, so dass ich unabhängig von der zeitlichen Länge des Signals, sondern nur in Abhängigkeit der Frequenz, das Ausgangssignal berechnen kann.
Geht das irgendwie? Bzw. kann ich mir aus H einen Filter basteln, damit ich diesen als direkt auf das Eingangssignal anwenden kann?
Was mir auch noch nicht ganz klar ist, in dem dspguide wird die obige Methode ja genauer beschrieben und darauf eingegangen, dass das Ausgangssignal nach der Rücktransformation die Länge N+M-1 hat, d.h. um M-1 Elemente länger ist als der Eingang, das ist bei mir aber gar nicht der Fall, wenn ich das Eingangssignal mit H (welches ich mit Nullen aufgefüllt habe) multipliziere.
Wäre dir denkbar, wenn du mir hier auf die Sprünge helfen könntest.
Ich habe nicht umsonst auf das n+m-1 hingewiesen...ein wenig mal selber darüber nachgrübeln würde sicher nicht schaden. Du must y auf die Länge n+m-1 entsprechend mit Nullen erweitern. In y ist der Ausschwingvorgang des Filters (m-1) nicht enhalten. Damit aber das Richtige rauskommt, muss erweitert werden. Wie groß du m wählst hängt letztendlich vom System ab. Wichtig ist, dass es innerhalb der Länge m eingeschwungen ist, falls ein IIR System. Ansonsten treten bei der Entfaltung Fehler auf. Mehr sage ich erstmal nicht dazu...denk mal drüber nach, warum man erweitern muss.
Die Impulsantwort h[m] bzw. Freq.antwort H[f] beschreiben genau dein System, dass aus einem Eingangssignal das Ausgangssignal erzeugt hat. Somit kannst du nun durch Faltung von h[m] und jedem beliebigen Testsignal ein Ausgangssignal erzeugen. Nichts anderes wird bei der Filterung gemacht...Faltung mit der Impulsantwort des Filters h[m].
abs(H[f]) + die richtige Skalierung ergibt dann das Betragsspektrum und zeigt dir somit das Frequenzverhalten deines Systems/Filters. Für welchen Freq.-bereich das gilt, hängt allein von der Abtastfreq. deines Eingangssignals ab
tfestimate
Gast
Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
Verfasst am: 24.05.2012, 00:11
Titel:
Das mit dem n+m-1 entsteht ja, soweit ich das richtig verstanden habe, durch die Faltung des Eingangssignals mit h(m).
Ich habe nun nochmal recht lange rumprobiert und weiter gelesen, allerdings komme ich immer noch nicht auf einen grünen Zweig...
Du schreibst, dass ich y auf die Länge n+m-1 mit Nullen erweitern muss.
Zur Bestimmung von H(f) benötige ich aber doch Vektoren gleicher Länge, um die Berechnung ohne Fehlermeldung durchführen zu können.
Wenn ich das Ergebnis von H(f) anschließend wie oben beschriebe überprüfe, erhalte ich auch das richtige Ergebnis. Wenn ich aber H(f) in h(m) umwandel und anschließend mit conv(in,h(m)) eine Faltung durchführe, erhalte ich y mit n+m-1 Elementen und wenn ich davon das Spektrum anschaue, ist es ein anderes als ich es durch Multiplikation von H(f) und fft(in) erhalte.
Heisst das also, dass ich schon bei der Bestimmung von H(f) einen Fehler mache, wenn ich x und y nicht auf n+m-1 Elemente erweiter?
Aus dem Skript hatte ich herausgelesen, dass ich, wenn H(f) und das Eingangssignal vorliegt und ich das Ausgangssignal bestimmen möchte, ich das Eingangssignal und H(f) entsprechend um Nullen erweitern muss, um der "circular convolution" zu entgehen. Damit erhalte ich aber ebenfalls keine sinnvollen Ergebnisse, vll aber auch, weil ich oben schon den Fehler mache.
Irgendwie scheine ich momentan auf dem Schlauch zu stehen und wäre daher für weitere Tipps dankbar, da ich nun auch nach zahlreichen Versuchen nicht wirklich weiter gekommen bin.
tfestimate
Gast
Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
Verfasst am: 24.05.2012, 11:19
Titel:
So nun habe ich endlich meinen Fehler gefunden
Ich habe ständig den Fehler gemacht, dass ich H(f) um Nullen erweitert habe, was natürlich absolut keinen Sinn macht. Natürlich muss ich h(m) um die entsprechende Anzahl von Nullen erweitern.
Anschließend kann ich y dann wieder auf die eigentliche Länge des Eingangssignals kürzen und erhalte die richtige Lösung.
Hier werden nun beide Signale auf die Länge n+m-1 gebracht. Denn wie du schon richtig erkannt hast, müssen die Signale die gleiche Länge aufweisen um multiplitziert oder dividiert werden zu können.
Nehmen wir an, du hast ein input und output Signal mit der Länge n = 1000. Dem Output Signal fehlt also somit die Länge m-1, denn sonst könnten durch die inv. Faltung nie wieder auf die 1000 Werte des Eingangssignals gerechnet werden. Wählt man nun m=100 und erweitert output nicht um 99 Werte, müsste ja das Eingangssignal ebenfalls 99 Messwerte weniger haben. Das Signal am Ende (m-1) wird in der Realität oft abgeschnitten, da es z.B. für die Filterung uninteressant ist. Die Funktion filter() in matlab macht es z.B. so. Tatsächlich liefert dir nur conv() diese m-1 Werte. Die m-1 Werte sind nach den 1000 input-Werten noch im Filter und müssten quasi durch nachschieben von Nullen aus dem Filter geholt werden. Damit aber die Faltung und inv. Faltung funktionieren, muss man unbedingt erweitern. Du hast hier auch noch ein zusätzliches Stichwort genannt...zyklische Faltung. Genau das ist das math. Problem.
tfestimate
Gast
Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
Verfasst am: 24.05.2012, 13:12
Titel:
DSP hat Folgendes geschrieben:
Und woher hast du h[m] wenn du eigentlich danach suchst?
h[m] hatte ich wie oben beschrieben durch ein gegebenes Eingangs- und Ausgangssignal gewonnen. So wie in deinem angegebenen Code, nur dass ich die Signale nicht erweitert hatte.
Obwohl ich diese Erweiterung bei der Berechnung von h[m] nicht angewendet hatte, sondern nur bei der Berechnung eines neuen Ausgangssignals h(m) und EIngangssignal auf die Länge m+n-1 gebracht hatte, habe ich die erwarteten Ergebnisse erhalten.
Kann dies damit zusammenhängen, dass mein Ein- und Ausgangssignal, mit welchem ich h[m] berechnet habe, schon von sich aus die gleiche Länge hatten?
Die Frage kurz zusammengefasst: Ist es für die Berechnung von h[m] nötig, das Eingangs- und Ausgangssignal zu erweitern, auch wenn diese Signale bereits gleich lang sind?
Lies dir bitte nochmal genau durch, was ich in meinem vorherigen Post geschrieben habe. Dort ist deine Frage nämlich schon beantwortet.
Außerdem kannst du mal folgendes Bsp. testen und dabei vor allem mal mit der Länge N spielen. Wenn du sie verkleinerst, enstehen Fehler bei der Rekonstruktion/Berechnung von h[m]. Ebenfalls spielt der Sprung als Eingangssignal für Sprungantwort eine Rolle....vgl. die 3 Zeilen für x_step.
Code:
clear;
%PT2 Strecke aufstellen und in z Bereich transformieren
%Streckenparameter
D1 = 0.5;
T0 = 0.0055;
Kst=1.00;
TS=1/4000; %Abtastzeit
N = 8192; % mal zum Test ändern
%zunächst im Bildbereich
s = tf('s');
%PT2 Strecke (= IIR Tiefpassfilter 2. Ordnung)
Gs = (Kst/(T0^2*s^2 + 2*D1*T0*s + 1));
G0 = Gs/(1+Gs);
% Diskretisierung
Gz = c2d(Gs,TS,'zoh');
% Zeitvektor
T = 0:TS:(N-1)*TS;
% Sprungantwort % Form des Eingangssprung spielt eine Rolle
%x_step = ones(N,1);
x_step = [zeros(100,1); ones(1000,1);zeros(N-1100,1)];
%x_step = [ones(1000,1);zeros(N-1000,1)];
% Impulsantwort durch Inverse Faltung von y_lsim/x_step
h_falt = inv_Faltung(y_step,x_step);
% Test: Rekonstruktion Eingangsignal durch inv. Faltung
x_rekonst = inv_Faltung(y_step,h_falt(1:N-1));
% Test: Ausgangssignal durch Faltung mit Impulsantwort
y_test1 = conv(x_step,h_lsim);
y_test2 = conv(x_step,h_falt(1:N-1));
figure(1);
hold on;
plot(T,y_step,'b*');
% m-1 abschneiden plot(T,y_test1(1:N),'g-.-');
plot(T,y_test2(1:N),'r--');
legend('step','Faltung mit orig. Impulsantwort','Faltung mit Rekonstruktion') title('Vergleich der Sprungantworten') hold off;
grid on;
% Vergleich Impulsantworten figure(2);
plot(T,h_lsim,'b');
hold on;
plot(T,h_falt(1:N),'r--');
legend('original','Rekonstruktion') title('Vergleich der Impulsantworten') hold off;
grid on;
% Vergleich des Betragsspekrums von realem Filter und Approximation
%w = (2*pi)./T;
[num, den] = tfdata(Gz,'v');
%h = tf(num,den,TS,'variable','z^-1');
[mag,f] = freqz(num,den,N,4000);
mag = abs(mag);
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.
RSS