Ãœbertragungsfunktion aus realen Signalen / Matlab / Octave ?

Hallo,
Markus Greim schrieb:
  Quoted Text Here  

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.
  Quoted Text Here  
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

Marcel Müller schrieb:
  Quoted Text Here  
D.h. möglicherweise bis in den Ural...
UP

Hallo Marcel,
  Quoted Text Here  
genau so...
  Quoted Text Here  
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.
  Quoted Text Here  
ok, das leuchtet ein, damit vermeidet man Fenster Effekte
  Quoted Text Here  
hast du dazu irgendeinen Literaturtipp ? Eventuell im Zusammenhang mit Matlab / Octave ?
  Quoted Text Here  
leuchtet ein
Grüße
Markus

Hallo,
Markus Greim schrieb:
  Quoted Text Here  
ja. Es hört sich exakt wie weißes Rauschen an.
  Quoted Text Here  
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.
  Quoted Text Here  
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.
  Quoted Text Here  
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

Hallo Marcel,
  Quoted Text Here  
dass wäre mal einen Versuch wert..
  Quoted Text Here  
soweit verstanden...
Die
  Quoted Text Here  
Latenzen hab ich eher nicht (denk ich ;-)) Die Ermittlung der Gruppenlaufzeit erfordert alledings ein
  Quoted Text Here  
gut
Danke für die vielen Tipps. Es scheint tatsächlich nichts "von der Stange" zu geben. Eigentlich unverständlich. Das sollte ja eigentlich eine elementare Ingenieursaufgabe sein, Das man da 60 Jahre nach Shannon noch bei Adam und Eva anfangen mus ist ziemlich frustrierend..
Grüße
Markus

Hallo Markus,
Markus Greim schrieb:
  Quoted Text Here  
als Alternative zu den anderen Vorschlägen hätte ich hier noch einen, der mir bei dem Thema Echokompensation aufgefallen ist, das mir gerade wieder über den Weg gelaufen ist. Dabei geht es um Freisprechen z.B. beim Telefon, um eine Rückkopplung des ankommenden Signals zu verhindern. Dazu wird ein adaptives Filter benötigt, das die Übertragungsfunktion des unerwünschten Pfades möglichst exakt nachbildet und die Rückkopplung unterdrückt:
http://en.wikipedia.org/wiki/Least_mean_squares_filter
Was dort als adaptives Filter bezeichnet wird, ist mir in Form eines FIR-Filters bekannt, dessen Koeffizienten durch eine Rekursion, eben die least mean square-Methode, berechnet wird. Wenn man ein lineares System ohne interference betrachtet, so ist diese Rekursion recht einfach durchzuführen.
Damit sich die Koeffizienten richtig einstellen, muß wohl das Eingangssignal x(n) alle Spektralanteile enthalten denke ich, also käme ein Sinussweep in Frage. Gemessen würde y(n) werden, und zwar _zeitgleich_ zu x(n). Wie gut und schnell das konvergiert und ob es dafür brauchbar ist, kann ich allerdings nicht sagen.
Mal angenommen, es konvergiert gut, dann ist ein Koeffizientensatz des (möglichst langen) FIR-Filters bekannt, der ja nichts weiter als die diskrete Impulsantwort des unbekannten Systems darstellt. Eigentlich müßte IMHO die Fouriertransformierte der FIR-Koeffizienten dann die komplexe Übertragungsfunktion sein. Die Auflösung im Frequenzbereich entspräche dann der (aus Symmetriegründen halben?) Länge des FIR-Filters.
Ob das Verfahren so brauchbar ist, könnte man in Matlab oder jeder beliebigen Programmiersprache testen, vielleicht erstmal an einem ganz einfachen bekannten System. Man erzeugt einfach den Datensatz für x(n), y(n), der der Messung entspräche und läßt damit die Rekursion durchlaufen, evtl. mehrmals. Wenn das Fehlersignal e(n)=0 wird, sieht es gut aus, denn dann muß das Filter ja stimmen. So sieht der Koeffizientenupdate aus:
e(k)^2=[y(k)-ye(k)]^2 Definition des Fehlerquadrats
x(k)=gemessenes Eingangssignal zum Zeitpunkt k y(k)= gemessenes Ausganggsignal des unbekannten Systems ye(k)=Ausgangssignal des adaptiven Filters e(k)=y(k)-ye(k)= Fehlersignal
FIR Koeffizientenupdate: Ci(k+1)=Ci(k)+u(k)*e(k)*x(k-i) 0<=u(k)<=2
Mit u(k) wird die Rückkopplung der Rekursion eingestellt, wäre experimentell zu ermitteln. Den Index i muß ich nachschauen, wo der anfängt und wie das FIR-Filter damit definiert ist.
mfg. Winfried

Hallo Winfried,
  Quoted Text Here  
ich denke das läuft so unter dem Oberbegriff Spektralschätzung.. Ich habe in den letzten Tagen noch etwas weiter recherchiert und bin auf eine schöne Veröffentlichung u. Software gestoßen die auch aus diesem Eck (Sprach Spektrogramme etc) kommt:
Time Frequency Toolbox for use with Matlab
von Auger, Flandrin, Goncalves u. Lemoine CNRS Frankreich u. Rice Univ. USA, 1995-1999
Das Tutorial findet man unter: http://tftb.nongnu.org/tutorial.pdf
und die Software und Homepage unter: http://tftb.nongnu.org /
Das ganze läuft unter dem Oberbegriff: "Analysis of non stationary signals using time frequency distributions"
weitere Stichworte sind: Wigner-Ville disribution, Short time fft, continous wavelet transformation...
Ich kann noch nicht endgültig sagen ob ich mit diesen Werkzeugen mein Problem lösen kann, es schaut aber bis jetzt sehr gut aus. Unabhängig davon ist das Tutorial unbedingt lesens- und und die Software (die auch unter Octave läuft) ausprobierenswert!
Danke soweit f. die vielen Anregungen.
Grüße
Markus Greim

Hallo Markus,
Markus Greim schrieb: [.....]
  Quoted Text Here  
der Index i müßte IMHO bei 0 anfangen und bis zum letzten FIR-Koeffizienten gehen. Mir ist noch eingefallen, daß das FIR-Filter ja eine Laufzeit hat, die im unbekannten System nicht vorhanden sein muß, wenn es keine Strecke mit ausschließlichem Echo ist. Dazu müßte man wohl eine Laufzeit als Ausgleich zum FIR in die Strecke legen. Diese bekannte Laufzeit kann man dann hinterher wieder rausrechnen.
Bei uns wurde diese Echokompensation mal in der Praxis getestet, dabei konvergierte es bei Rauschen gut und bei Sprache weniger gut. Dabei ging es aber wirklich nur um Echos, ob das Ganze auch bei beliebigen Frequenzgängen konvergiert, kann ich nicht sagen, wäre zu testen. Von der Plausibilität her scheint mir der Laufzeitausgleich zwingend zu sein, da FIR-Filter immer laufzeitbehaftet sind im Vergleich zu IIR-Filtern. Der Vorteil dieses Verfahrens, wenn es funktioniert, wäre die extreme Einfachheit.
  Quoted Text Here  
Ich habe es mir mal runtergeladen, vielen Dank. Mit der Thematik Wavelet habe ich mich noch nicht beschäftigt, sieht sehr theoretisch aus. Eigentlich müßten die Leute aus der Theorie der Regelungstechnik was dazu wissen, denn das müßte ja deren Job sein, also die Erkennung eines Systems.
  Quoted Text Here  
Viel Erfolg! Vielleicht spiele ich auch mal damit herum, wenn ich Lust und Zeit habe.
mfg. Winfried