Übertragungsfunktion aus realen Signalen / Matlab / Octave ?

Hallo,

irgendwann holt einen die Systemtheorie wieder ein....

Mein Problem: gegeben sei:

- ein Einganssignal (E(t_i)) (abgetastet im Zeitbereich) , Typ sinus sweep, konstante Amplitude, Frequenz 0..x Hz.

- ein Ausgangssignal als Antwort auf das Eingangssignal. Das System sei zeitinvariant. (A(t_i))

gesucht: Real- und Imaginärteil der Übertragungsfunktion S(f)

was ich weiß:

S(f) = Wirkung / Ursache = A(f) / E(f)

was ich nicht weiß : wie man S(f) oder S(f_i) also diskret numerisch am geschicktesten berechnet.

Ich hab mal naiv in Octave (Matlab clone)

fft(Ai)./fft(Ei)

gerechnet, aber ich glaub nicht das das stimmt.

Es geht hier nicht um eine Theorie oder Hausaufgabe sondern um praktische Messtechnik. Diverse Fundstellen im web scheinen alle immer nur die Theorie voneinander abzuschreiben. Als Literatur hab ich mal R.Unbehauen, Marko und Stearns konsultiert. Das sind aber alles eher Bücher aus meiner eher vordigitalen (Studien)-Zeit (ich bin nicht mehr der Jüngste :-) )

Irgendwelche Hinweise oder wenigsten treffende Literaturhinweise ??

Danke und Grüße

Markus Greim

Reply to
Markus Greim
Loading thread data ...

Hallo, Markus,

Du (greim) meintest am 17.07.08:

Also musst Du zuerst (z.B. mit Hilfe der Fourier-Transformation) aus dem Datenhaufen so etwas wie ein Bündel von sinusförmigen Schwingungen machen.

Wenn ich die Kurzbeschreibung richtig verstanden habe: das Eingangssignal ist (pro Messung) monofrequent und stabil ("reiner Sinus")?

Wie wäre es, wenn Du zuerst nur die Fouriertransformation machst und Dir die dort angebotenen Ergebnisse anschaust?

Da solltest Du zuerst prüfen, ob das Bündel von Frequenzen "sinnvoll" aussieht und welche Frequenzen Du getrost weglassen kannst (nennt sich "Tiefpass").

Viele Gruesse! Helmut

Reply to
Helmut Hullen

"Markus Greim" schrieb im Newsbeitrag news:g5n7cr$nqd$02$ snipped-for-privacy@news.t-online.com...

E(s) = 1 Einheitssprung S(s) = A(s) Sprungantwort (Übertragungsfunktion)

Siehe Approximation Seite 5

formatting link

Reply to
JCH

JCH schrieb:

jou, so weit war ich schon. Aber bitte siehe unter "gegeben": Ich kann auf das (mechanische) System keinen Sprung aufbringen, das Sinus Signal ist zwingend!

Markus

Reply to
Markus Greim

Markus Greim schrieb:

Das ist die klassische Verwendung eines Vektor-Netzwerkanalysators.

Du bekommst für jede Frequenz Amplitude und Phase, was beliebig in Real- und Imaginärteil umgerechnet werden kann.

Unter der Voraussetzung, dass die Sweepgeschwindigkeit angemessen ist, kannst Du das Eingangssignal auf Nullstellen untersuchen, somit eine ganze Periode auffinden. Das Ergebnis ist bei einer linearen Übertragungsfunktion ebenfalls eine Sinuswelle, deren Amplitude und Phase zu in Bezug auf das Eingangssignal setzen kannst. Eine FFT ist hier eigentlich nicht angebracht oder besser gesagt erhältst Du ohnehin unter den o.a. Voraussetzungen nur den Gleichspannungsanteil und die Grundwelle.

Gruss Udo

Reply to
Udo Piechottka

Moin,

Markus Greim schrub:

Du musst dir mal ansehen, wie Octave das Ergebnis der FFT im Ergebnisvektor überhaupt ablegt. Das muss nicht unbedingt auf naheliegende Weise passieren. An sonsten sollte es richtig sein. Allerdings: Der Weg mit der FFT geht davon aus, dass es sich um ein zyklisches Signal handelt. Durch diese Annahme können Schmutzeffekte entstehen. Am besten: Starte die Messung bei ausgeschalteter Maschine und beende sie auch bei ausgeschalteter Maschine. Zwischendurch einschalten und den Sweep machen. Also kurz gesagt eine _lange_ Messreihe aufnehmen im Vergleich zu den zu beobachtenden Frequenzen.

Was natürlich noch die Sache stören kann: Dein System könnte nichtlinear sein. In dem Fall wirst du sicher sehr komische Ergebnisse bekommen. Praktisch ist jedes mechanische System nichtlinear. Deswegen solltest du das etwas in Hinterkopf haben, wenn du das Ergebnis betrachtest: Tritt z.B. bei 50Hz ein starker Durchgriff (großes Ausgangssignal) auf und ebenso bei

100Hz, 150Hz,... dann haben die Werte bei 100Hz,... eventuell nicht viel zu bedeuten, es sind nur Artefakte aus der Nichtlinearität bei der 50Hz-Schwingung.

CU Rollo

Reply to
Roland Damm

Hallo,

Markus Greim schrieb:

doch, das mache ich andauernd so (gesetzt das / eine Komponentenweise, Komplexe Division ist). Bei Sweep muss man aber ein FFT über die Gesamte Aufnahme machen. Das kostet mächtig Speicher und CPU. Ferner wird man pro Frequenzkanal das SNR genau betrachten müssen. Nur da wo fft(Ei) signifikant Amplitude hat, kommt auch etwas vernünftiges raus. Und da hat mna bei einem Sweep das Problem, dass natürlich jede Frequenz nur kurz dran ist und die meiste Zeit sich also nur Rauschen auf diesem Frequenzkanal sammelt. Ich nehme dafür für diese Messart entweder MLS-Signale oder definierte kümstliche Rauschsignale die durch ifft erzeugt werden.

Sweep ist nicht so der Brüller. Wenn man den Sweep in der Frequenz quantisieren könnte, wäre das besser. Also N diskrete Frequenzen, die allesamt ein Vielfaches eines geeignet zu wählenden Messintervalls sind das einer Zweierpotenz an Samples entsprechen sollte. Dann kann man folgendes machen:

Generator Messstatus

------------ ----------------------- Settle Delay Frequenz 1 ----------------------- Sample für Frequenz 1

------------------------------------- Settle Delay Frequenz 2 ----------------------- Sample für Frequenz 2

------------------------------------- ...

Die Sampledauer muss immer ein vielfaches der Periodizität der zu messenden Frequenz sein. Dann berechnet man für jede Frequenz einzeln fft(Ea)/fft(Ei) nimmt aber nur den Koeffizienten des Ergebnisses, der der aktuellen Anregungsfrequenz exakt entspricht. Damit hat man so eine Art Lock-In gemacht und erreicht selbst unter räudigen Bedingungen exzellente SNRs (z.B. >60dB mit einem 8-Bit ADC in störungsreicher Umgebung). Diese Qualität hat man nicht, wenn man die Messmethode nicht mit dem Wissen anreichert, welche Frequenz gerade gemessen wird. Mathematisch gesehen kann man sich die FFTs natürlich sparen, wenn man nur einen Koeffizienten braucht und einfach eine PCA mit einem Quadratursignal der Anregungsfrequenz machen. Das ist extrem schnell, da O(n) in CPU und O(1) im Speicher.

Theoretisch bekommt man fast genausogute Ergebnisse, wenn man einen kontinuierlichen Sweep geschickt auswertet. Aber das ist /erheblich/ mehr Aufwand und birgt Risiken, wie beispielsweise Schatten hinter Polen mit hoher Güte.

Marcel

Reply to
Marcel Müller

Marcel Müller schrieb:

D.h. möglicherweise bis in den Ural...

UP

Reply to
Udo Piechottka

Hallo,

Das wäre natürlich ideal, aber Hardware kann ich hier nicht einsetzen. Ich muss die beiden Datensätze E_i und A_I softwaremäßig auf einem Linux Rechner verarbeiten

Grüße

Markus

Reply to
Markus Greim

Hallo Marcel,

genau so...

MLS ist mir neu. Hab gerade etwas danach gegoogelt Wie "hört" sich denn so ein MLS Signal an? Ist das ein spezielles rauschen ? Die Frage ist ob man ein mechanisches System damit anregen kann. Mein Aktor kann nur begrenzt Energie / Zeiteinheit umsetzen.

ok, das leuchtet ein, damit vermeidet man Fenster Effekte

hast du dazu irgendeinen Literaturtipp ? Eventuell im Zusammenhang mit Matlab / Octave ?

leuchtet ein

Grüße

Markus

Reply to
Markus Greim

Hallo, Markus,

Du (greim) meintest am 18.07.08:

In aller Unkenntnis der speziellen Daten: FFT sollte eine Summe von Daten liefern (für diverse Frequenzen), und Summen bei Brüchen sind Teufelszeug.

Da würde ich gern nachfragen - bei Deiner Vermutung ergibt sich ja kein kontinuierliches Spektrum, sondern ein diskretes (okay - die Frequenzen sollten ganzzahlig sein).

Ist "sweep" das, was früher mal "wobbeln" war? Alle Frequenzen von 1 Hz bis 1 GHz kontinuierlich anlegen?

Viele Gruesse! Helmut

Reply to
Helmut Hullen

Markus Greim schrieb:

Daran hindert Dich niemand. Aber eine FFT über einen beliebigen Sweep laufen zu lassen, ist ungünstig, vor allem wenn Du doch schon weisst, dass es ein Sinussignal ist und Du die Frequenz ermitteln kannst.

Ich meinte schon, dass Du das numerisch lösen wirst. Du kannst beispielsweise die den Momentanwert der Frequenz aus deinem Datensatzermitteln und hast somit eine volle Periode, die Du analysieren kannst. Du hast somit den stationären Zustand bei einer Frequenz und kannst Amplitude und Phase bestimmen. Per Euler lässt sich das in Re und Im umrechnen.. fertig.

Und das machst Du jetzt für jede Periode oder mit entsprechender Reduktion nur über jede n-te Periode und bekommst somit den frequenzabhängigen Verlauf der gesuchten Grössen, also den Frequenzgang.

-Udo

Reply to
Udo Piechottka

Hallo Helmut,

genau, hier liegt glaube ich das Problem, wie dividiert man diese beiden Summen numerisch am geschicktesten ????

ja, ja ich hab nur gedacht so altmodische Begriffe verwendet man nicht mehr ;-) Da wir es mit einem mechanischen System zu tun haben ist der Frequenzbereich eher 0..20Hz

Grüße

Markus

Reply to
Markus Greim

Hallo, Markus,

Du (greim) meintest am 18.07.08:

Zerlege es in Teilschritte:

Anrege-Frequenz 5 Hz, Ergebnisdaten Anrege-Frequenz 10 Hz, Ergebnisdaten usw.

Dann ist einer der beiden Teile des Bruches keine Summe mehr.

Also diskretisieren. Und dann schaust Du Dir die Ergebnisse der einzelnen Anrege-Frequenzen an und versuchst sie zu vereinfachen.

Viele Gruesse! Helmut

Reply to
Helmut Hullen

Helmut Hullen schrieb:

Das will er offenbar nicht, sucht möglichgerweise eine fertige Funktion in der Bibliothek... Habe ihm auch schon vorgeschlagen, periodenweise und dn sinnvollen Abschnitten diskret zu rechnen...na denn

Gruss Udo

Reply to
Udo Piechottka

Hallo Udo, Hallo Helmut,

darauf wird es wohl hinauslaufen. Ich dachte nur es gäbe irgendeine elegantere Technik

Grüße

Markus

Reply to
Markus Greim

Hallo,

Markus Greim schrieb:

ja. Es hört sich exakt wie weißes Rauschen an.

Die Frage ist, welche Variable dem Signal folgt. Wenn es die Beschleunigung ist, geht das natürlich. Bei Geschwindigkeit oder gar Ort ist es Unsinn. MLS-Signale sind halt einfach zu erzeugen, haben einen Crest-Faktor von 1 und damit eine bei gegebener Maximalamplitude optimale Energiedichte. Mehr nicht. Wenn man Ortskurve eines mechanischen Systems aufprägen möchte, tut man im allgemeinen sowieso gut daran, kein weißes Rauschen zu verwenden, sondern Rosa Rauschen oder ähnliches.

In anderen Fällen nehme ich gerne periodisches Rauschen. Also im einfachsten Fall ein stück Rauschen immer wiederholen. Bei Weißem Rauschen funktioniert das einfach so, und man bekommt statt einem kontinuierlichen Spektrum eben ein diskretes. Im allgemeinen Fall fährt man mit folgender Generatorfunktion für bandbreitenbegrenztes Rauschen gut:

Ei(t) = ifft( if (f element ]fmin, fmax]) then (f^kappa * exp(2*pi*i*random()) else 0 )

Das Ergebnis normalisiert man dann einfach auf ein für die Anwendung sinnvolles Eimax.

kappa=0 und fmin=0 und fmax=fNyquist gibt weißes Rauschen ohne Bandbreitenlimit. kappa=-1 gibt Rosa Rauschen. fmin muss dann natürlich sinvoll gesetzt werden.

Je nach Anwendungsfall kann man kappa variieren. Sinvolle Werte liegen typischerweise zwischen -1,5 und 0,5 in 0,5er Schritten.

Solange man keine systembedingte Integration oder Differentiation in der zu messenden Übertragungsfunktion hat und man eine lineare Frequenzskala im Ergebnis haben will, tut man gut daran, kappa auf 0 zu lassen. Für eine logarithmische Frequenzskala setzt man kappa um eins herab. Dann muss man allerdings in der Auswertung auch benachbarte Frequenzkanäle im Ergebnis der FFTs zusammenfassen, so dass man effektiv Frequenzkanäle ähnlicher logarithmischer Breite bekommt. Die Zusammenfassung muss sinnvollerweise durch Mittelwertbildung im Endergebnis fft(Ea)/fft(Ei) und in Amplitude und Gruppenlaufzeit erfolgen. Im Endergebnis deshalb, weil die Phasen benachbarter Kanäle im Anregungssignal zufällig, also unkorreliert sind. Über die Gruppenlaufzeit deshalb, weil die Phase der Übertragungsfunktion üblicherweise schnell veränderlich sein kann, wenn man es mit Latenzen zu tun hat. Die Ermittlung der Gruppenlaufzeit erfordert alledings ein Phase-Unrolling. Wenn die Frequenzstützpunkte dicht genug liegen, kann man einfach die in Polarkoordinaten (Ammplitude, Phase) kürzere Entfernung nehmen. Größere Latenzen sollte man zuvor per Multiplikation mit exp(2*pi*i*f*linearshift) kompensieren. Theoretisch kann man linearshift per Least-Square-FIT bestimmen, aber das ist wegen der Zahlreichen lokalen Minima ziemlich instabil. Schneller geht es mit der Hand und typischerweise ist es sowieso eine Konstante des Messaufbaus. Genau stimmen muss der Wert ohnehin nicht, da es keinen unmittelbaren Einfluss auf das Ergebnis gibt.

Ja, die Fenstereffekte werden genau dann vermieden, wenn das Anregungssignal in einer festen (aber unbekannten) Phasenbeziehung zu den aufgenommenen Ei- und Ea-Signalen steht und die FFTs für Testsignalgenerierung und Analyse gleich lang sind. Dazu müssen die Sampleratengeneratoren der Ein- und Ausgänge aus demselben Quarz synchronisiert werden. Das Kriterium erfüllen glücklicherweise die allermeisten Messkarten und selbst billigste Soundkarten. Ich würe Ei auch immer nochmal mit aufzeichen und ausschließlich dieses Signal zur Auswertung verwenden. Damit sind jegliche Latenz- und Dreckeffekte der Messkarten und auch des Aktors weitestgehend kompensiert (typisch -50dB), und man ist nur noch sensitiv auf die Differenzen zwischen den beiden (hoffentlich baugleichen) Aufnehmern.

Nein, sorry. Ich habe mir sowohl das Verfahren als auch die Software selber zusammengepuzzelt. Die Auswertung mache ich einfach mit zwei kleinen C-Programmen. Eines liefert das Testsignal an stdout, was ich dann per Pipe an das PCM-Sounddevice schicke. Das andere bekommt Ea und Ei von PCM-Device (Line-In, Stereo) an stdin und führt die notwendigen Berechnungen durch, schreibt nach jedem FFT-Zyklus (bei mir zwischen

65536 und 262144 Samples) ein Ergebnisfile mit der Übertragungsfunktion und gibt über stdout das Visualisierungskommando an gnuplot. Optional mittle ich über alle bisher gemessenen Perioden. Die Mittelung findet allerdings auf den rohen Samples statt, da alle Nutzsignalanteile notwendigerweise periodisch sind.

Im (stufenweisen) Sweepmodus synchronisiere ich einmal am Anfang mit einem Testsignal mit Phasensprung und zähle ab da in beiden Programmen parallel die Samples um die verschiedenen Frequenzschritte hinter sich zu bringen. Die Messung kann dann bis zu einer Dreiviertelstunde alleine laufen. Sweep ist hochpräzise aber langsam. Mit Rauschen hat man meist schon nach einigen Sekunden die ersten brauchbaren Ergebnisse. Hüten sollte man sich vor der Rauschmessung, wenn man mit Nichtlinearitäten rechnet oder ohnehin Probleme mit dem SNR hat.

Ich mache auf die Weise unter anderem auch Audiomessungen von Lautsprechern. Im Sweepmodus kann man problemlos im Garten mit ordentlicher Geräuschkulisse messen - upfire und schräg neben der Hausecke versteht sich, damit nichts zurück kommt. Die Ergebnisse passen exakt über die Messkurven professioneller Labors mit schalltoten Räumen und sind meist sogar weniger verrauscht und reicher an Stützpunkten. Aber genausogut kann man auf die Tour andere Netzwerkanalysen (Zweipol und Vierpol) machen oder ein präzises LCR-Meter (bei moderaten Frequenzen) basteln. Mit einem DSO kann man auch in gehobene Frequenzbereiche vordringen. Manche davon beherrschen sogar die notwendigen Auswertungsschritte FFT Amplitude, FFT Phase, Amplitudenquotient und Phasendifferenz on the fly. Allerdings ist es nicht trivial, ein synchrones Testsignal zu bekommen, das selbst professionelle DSOs meines Wissens den Wandlertakt nicht rausrücken.

Marcel

Reply to
Marcel Müller

Markus Greim schrieb:

So etwas selber zu programmieren ist anstruchsvoll, erfreut aber das Gemüt wenn es dann funktioniert und das auch noch elegant.

Gruss Udo

Reply to
Udo Piechottka

"Markus Greim" schrieb im Newsbeitrag news:g5o0ak$4ln$01$ snipped-for-privacy@news.t-online.com...

O.g. 'Beispiel' zeigt eine Approximation für Dfgl 5. Ordnung mit starker Dämpfung.

Ein mechanisches Beispiel könnte z.B. ein Pendel mit Dämpfung sein:

T^2*s^2 + 2*d*T*s + 1 = C

d
Reply to
JCH

Hallo JCH ,

das ist mir klar, wenn ich die Üfkt. denn hätte!

wenn ich deine Ausführungen richtig interpretiere empfiehlst du mir für die Strecke ein Modell aufzustellen, und dann die Parameter dieses Modelles zu approximieren ??

Grüße Markus

Reply to
Markus Greim

PolyTech Forum website is not affiliated with any of the manufacturers or service providers discussed here. All logos and trade names are the property of their respective owners.