Übertragungsfunktion aus diskreten Werten

https://www.mikrocontroller.net/articles/Digitalfilter_mit_ATmega#Frequenzgang_ma.C3.9Fgeschneidert
speziell der den Artikel zu FDLS:
http://www2.units.it/ramponi/teaching/DSP/materiale/Ch5%28x%293e.pdf#page=29
"the FDLS algorithm produces a transfer function that approximates an 
arbitrary frequency response"
also aus Messwerten wird eine Übertragungsfunktion (ein IIR-Filter 
wählbarer Stufenzahl) angenähert.
Nein, das war die falsche Methode da müssen ja unterschiedliche 
Frequenzen gemessen werden also im Frequenzbereich, dann ist es 
vielleicht der "prony"-Algorithmus

Du erhälst mit der Division der Fouriertransformierten des 
Ausgangssignals durch die Fouriertransformierte des Eingangssignals die 
Fouriertransformierte der Impulsantwort. Das hat bei mir, als ich das 
damals gebraucht habe, ganz gut funktioniert.

Eine weitere Möglichkeit wäre die Verwendung eines adaptiven Filters mit 
LMS, NLMS oder RLS Algorithmus.

Das alles geht aber nur, wenn dein System linear ist!

Wenn man diskrete Werte hat und eine Funktion davon moechte, muss man 
ein rationales Polynom fitten.

Das Einfachste, aber Aufwendigste ist N Gleichungen mit N unbekannten. 
Es kann sein, dass das System einen zu hohen Grad aufweist und die 
Koeffizienten nicht stabil sind. Das System ist dann schlecht 
konditioniert.

Dann macht man ein einfacheres Modell, das nun aber ueberbestimmt ist. N 
Gleichungen mit M Unbekannten, mit N>M. Dann muss man einen Least Square 
Fit machen.

Ich wuerde das zweite verwenden, auch wenn es mathematisch aufwendiger 
ist.

Noch einfacher ist wenn man schon eine bessere Vorstellung hat, dann 
kann man die Funktion graphisch fitten. Eine Programmieruebung.

Hallöchen, seid ihr teilweise kirre?

Sorry Leute, aber man differenziert keine Sprungantwort, um daraus ne 
Impulsantwort zu bekommen, viel zu rauschempfindlich (Differentiation 
ist ein super Hochpass).
Ebenso soll hier kein Polynom gefittet werden, das ist ein dynamisches 
System, nix statisches.

Das Zauberwort heißt Systemidentifikation, wurde oben (zwischen den 
absurden Antworten) schon genannt.

Dazu passt Kapitel 8 "Methode der kleinsten Quadrate für dynamische 
Prozesse" aus dem Isermann, "Identifikation dynamischer Systeme 1". 
Davon gibts noch diverse Modifikationen (rekursiv, gewichtet, totale ls 
usw usw).

Du möchtest also eine Übertragungsfunktion eines zeitdiskreten Modells 
haben?

mit dem Vektor der Eingangssignale u(z), dem Vektor der gemessenen 
Ausgangssignale y(z), der Ordnung m und N der Anzahl der Samples (im 
Buch kommt noch eine Totzeit dazu, lasse ich der Einfachheit halber 
jetzt weg)

Sind einige Voraussetzungen erfüllt:
- Ordnung m und Totzeit d bekannt oder zumindest grob bekannt (immer 
schön niedrig anfangen...)
- Eingangssignal muss Messbar sein (ist erfüllt)
- Störsignal ist mittelwertfrei und stationär (z.B. weißes Rauschen 
durch die Messungen). Kp ob das erfüllt wird bei dir, geht aber häufig 
auch ohne ;-)

Dann kannst du deine Parameter ai und bi schätzen, wenn du

 rechnest.
Mit

,

, du stapelst also einfach alle deine Messungen übereinander. Dann 
brauchst du noch die Regressormatrix

, stapelst also auch deine Eingangs und Messgrößen in eine Matrix.
Dann berechnest du wie oben angegeben die Pseudoinverse und hast deine 
Parameter ai und bi.


Weitere Stichworte sind Parameteridentifikation, output error model, 
least squares, armax, system identification. Das ist alles in der System 
Identification Toolbox in Matlab eingebaut (und noch viel viel viel 
mehr).

Schöne Grüße,
Jan

Hallöchen Jan,
anzumerken bei deinem Ansatz wäre noch, dass dies nur für 
nicht-rekursive Systeme funktioniert. Der Grund ist der, dass im Vektor 
y das gemessene Ausgangssignal ist, die Pseudoinverse minimiert also den 
Fehler, wenn an deinem System genau dieses Ausgangssignal raus kommt. 
Leider ist das Ausgangssignal aber bei einem rekursiven System eine 
Funktion vom Eingangssignal und von sich selbst (zeitverschoben). Ein 
IIR-System kannst du so nicht identifizieren.



Jan K. schrieb:
> Hallöchen, seid ihr teilweise kirre?
>
>
> Das Zauberwort heißt Systemidentifikation

Uhoh, da sind leider 1, 2 Fehler drin sorry.
u(m) in der vorletzten Spalte ist nicht richtig, da muss natürlich 
u(m-1) hin. Auch y ist falsch wie du richtig erkannt hast, y muss bei 
y(m) starten, nicht bei y(0).
Die Gesamtanzahl deiner Samples ist N+m.

Essenziell für dein Problem ist übrigens, dass die Werte vorher mit 
einer konstanten Abtastrate resamplet werden (dein diff(t) ist nicht 
konstant!), sonst funktioniert das ganze nicht.

Ich bekomme übrigens

für eine first order approximation heraus. Scheint nicht zu reichen, 
siehe Bild #3 ;-)

@Tobias
Doch, man kann auch rekursive Systeme mit obigem Ansatz identifizieren:

Das ist ein reines IIR System 3. Position der Pole und Nullstellen kann 
man im ersten Bild sehen.

1

[y, tt] = stepz(1,[1 -.95 .001 .1],100); % daten generieren

2

u = ones(size(y)); % eingangssignal

3

m = 1; % approximation 1. Ordnung

4

psi = [cell2mat(arrayfun(@(i)-y(m-i+1:end-i),1:m,'unif',false)) cell2mat(arrayfun(@(i)u(m-i+1:end-i),1:m,'unif',false))]; % regressormatrix

5

p = psi\y(m+1:end); % parameter

6

7

figure,

8

stem(y)

9

hold all

10

stepz([0 p(floor(length(p)/2)+1:end)'],[1 p(1:floor(length(p)/2))'])

11

legend('original','approximated')

First order approximation resultiert also in einem Modell

Beachte, dass im Nenner der erste Koeffizient auf 1 normiert ist (wird 
nicht mit identifiziert) und dieser im Zähler == 0 ist.

Du kannst auch die Modellordnung höher wählen, nur dann ist das 
Gleichungssystem nicht mehr gut konditioniert. Funktionieren tut es 
trotzdem, probier es mal mit m=2 oder m=3.

@asma
wie bist du auf die positiven Exponenten gekommen? Wenn ich das umrechne 
bleibt nur ein konstanter gain im Zähler übrig, das kann diese Methode 
eigentlich nicht berechnen.

@Jan

man muss vielleicht präzisieren:
man kann ein IIR-System so identifizieren und bekommt dann eine 
Übertragungsfunktion, welche den Gleichungsfehler (so hiess das glaube 
ich) minimiert. In der Regel möchte man aber den Ausgangsfehler 
minimiert haben, und das geht so nicht!

Grund:
dein Vektor y enthält die gemessenen Werte des Ausgangssignals. Wenn das 
ein IIR-System sein soll, dann sind das aber auch gleichzeitig 
Funktionen der Koeffizienten a. Da du die a aber nicht kennst, .... du 
siehst das Problem.

Ich bin mir nicht sicher, ob "Oh Doch (jetztnicht)" das meinte, oder ob 
er einfach nur ein statisches Polynom fitten wollte.
Das Modell ist zwar polynomial, aber eben im Kontext einer 
Übertragungsfunktion.

Woher hast du deine Quelle? :-)

Da ist die Namensgebung sehr blöd von mir und beruht auf Faulheit :-/

Wenn y alle deine m+N gemessenen Werte sind, so sei Y ein Vektor mit 
allen Messungen, bei y(m-1) beginnend. Y ist dann der Vektor auf der 
rechten Seite des LGS oben.

Versuche das Matlab Skript oben nachzuvollziehen.
Bei p = psi\y(m+1:end); wird nicht der gesamte Vektor y, sondern eben ab 
dem m+1-ten Element genommen.

Wichtig: Die Schreibweise in der Matrix oben geht von der Indizierung 
von 0 aus, y(0) ist also dein erstes gemessene Ausgangselement zum 
Zeitpunkt k=0. In Matlab muss alles um +1 verschoben werden, da das 
erste Element am Index 1 ist.

Wat bist du denn für einer? Danke, dass du dich um mein Privatwohl 
kümmerst.

Nein, man macht natürlich kein Foto, schonmal was von Copyright gehört? 
Außerdem ist es in dem Buch noch etwas komplizierter dargestellt.

Die Matlab Version funktioniert und sollte mit der Matrix oben 
übereinstimmen.

Wie das LGS gelöst wird steht oben schon. Matlab bildet beim "\" bzw 
"mldivide" automatisch die Pseudoinverse [bzw berechnet die implizit 
über eine QR Zerlegung], aber man kann das natürlich auch zu Fuß 
ausrechnen. Wie gesagt, wie die Pseudoinverse definiert ist steht in 
meinem ersten Post. Wie man Matrizen transponiert und invertiert wissen 
hoffentlich alle, wenn nicht ist der Gauß-Jordan Algo eine erste 
Anlaufstelle.

Bin mir bei dir nicht sicher, ob du trollst oder einfach irgendwo was 
gegoogelt hast.