mikrocontroller.net

Forum: Digitale Signalverarbeitung / DSP Übertragungsfunktion aus diskreten Werten


Autor: derneue (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Ich habe diskrete Werte (Zeit und dazugehöriger Wert) die aus einem 
Einheitssprung entstanden sind. Wie bekomme ich daraus nun die 
Impulsantwort/Sprungantwort oder die Übertragungsfunktion?

Geht da was mit Laplace? Oder mit Hilfe von Matlab? Kenne mich mit 
Regelungstechnik leider nicht aus.

: Verschoben durch Moderator
Autor: aSma>> (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Wenn du die Werte schon hast:
-Excel
-Gnuplot
-Matlab
such dir was raus.

Autor: derneue (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Okay, aber wie genau was muss ich machen?

Autor: aSma>> (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Kommen jetzt alle Leute hier ausm Kindergarten?!
-Googeln. Soll man das für dich machen?
Vorallem du sagst nicht einmal was du überhaupt für ein Programm zur 
Vergügung hast. Viel Spaß.

Autor: derneue (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Habe zugriff auf Matlab und MS Excel. Gegoogelt habe ich bereits aber 
nichts treffendes gefunden deshalb frage ich hier.

Autor: Christoph Kessler (db1uq) (christoph_kessler)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
https://www.mikrocontroller.net/articles/Digitalfi...
speziell der den Artikel zu FDLS:
http://www2.units.it/ramponi/teaching/DSP/material...
"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

: Bearbeitet durch User
Autor: aSma>> (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
derneue schrieb:
> Ich habe diskrete Werte (Zeit und dazugehöriger Wert) die aus einem
> Einheitssprung entstanden sind. Wie bekomme ich daraus nun die
> Impulsantwort/Sprungantwort oder die Übertragungsfunktion?

Sry, ich habe dich falsch verstanden. Du hättes mal die Folge aufmalen 
sollen. Letztendlich musst du nur den Quotienten bilden:

w(z)=1z^0 + 1z^-1 + 2^-2 (Eingangsprung)
x(z)=4z^0 + 3z^-1 + 2^-2 (z.B. gegebene Folge)
G(z)=x/w jetzt nur noch mit den größten Exponent erweitern.
Hier mit z^2/z^2.

Im Simulink braucht man dies nichteinmal zu machen. Einfach einen zwei 
Blöcke mit Eingangsprung und Transfer Function erstellen + scope.

In Matlab ist dieser Link nützlich:
https://de.mathworks.com/help/symbolic/compute-z-t...
Da muss man erstmal die Übtragung in den Zeitbereich transformieren.

mfg

Autor: derneue (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
aSma>> schrieb:
> w(z)=1z^0 + 1z^-1 + 2^-2 (Eingangsprung)
> x(z)=4z^0 + 3z^-1 + 2^-2 (z.B. gegebene Folge)
> G(z)=x/w jetzt nur noch mit den größten Exponent erweitern.
> Hier mit z^2/z^2.
>
> Im Simulink braucht man dies nichteinmal zu machen. Einfach einen zwei
> Blöcke mit Eingangsprung und Transfer Function erstellen + scope.

Ja abe wie kommst du auf die Gleichungen, die sehen bei mir bestimmt 
anders aus? Ich weiss leider immernoch nicht wie ich vorzugehen habe. 
Simulink ist auch adrauf bei Matlab bei mir wie könnte ich nun 
vorangehen. Ich kann meine Tabelle mit den Werten in eine cell variable 
oder struct variable laden. Und dann?

Autor: Der Andere (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
derneue schrieb:
> Geht da was mit Laplace?

Nennt sich im diskreten Z-Transformation. Ist aber schon fast 30Jahre 
her bei mir, ich müsste genauso nachlesen wie du. Also Suche nach 
Z-Transformation

Autor: aSma>> (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
derneue schrieb:
> Ja abe wie kommst du auf die Gleichungen, die sehen bei mir bestimmt
> anders aus?

Mit Verlaub, ja. Stell deine Folge hier rein, dann sehen wir weiter.

> Ich weiss leider immernoch nicht wie ich vorzugehen habe.
Habe ich dir schon zuvor an einen Beispiel gezeigt.

Autor: derneue (Gast)
Datum:
Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Links die spalte ist die Zeit rechts die dazugehörigen Werte.

Autor: aSma>> (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Ich nehme alles wieder zurück. Du bist doch unfähig!!!

Speicher die deine datei in den current folder als: werte.dat
dann tippst du das rein:
>> load werte.dat
>> plot(werte)

Autor: derneue (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Und was soll das? Ich weiss wie die Kurve aussieht, ich weiss auch wie 
ich sie plotten kann. Ich brauche doch die Übertragungsfunktion aus 
diesen Daten.

Autor: aSma>> (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hilf dir das weiter ab S.8
http://staff.ltam.lu/feljc/school/asser/2_Regelstrecken.pdf

Man kann eine Pt2 Strecke als Pt1-Strecke mit einer Totzeit nachbilder 
oder durch einen allpass glied aus Pt1 strecken.

Oder willst du das digital bearbeiten können?

Autor: Alexander (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Ich glaube, dass es ihm darum geht, die Pol- und Nullstellen des Systems 
anhand einer aufgenommenen Folge zu ermitteln.

Das Stichwort heißt hier: Systemidentifikation.

Ich habe mir die Messungen nicht angesehen, wenn es allerdings eine 
zeitliche Messung ist (Sprungantwort etc.), bietet sich z.B. die Methode 
der kleinsten Fehlerquadrate (Least Square) an.

Gruß,

Autor: aSma>> (Gast)
Datum:
Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Hier das ist eine PT1 Strecke. Die geringe Verzugszeit ist eigentlich 
vernachlässigbar.
clear all;
clc;

load data.dat

s = tf('s');
K = 100;
Ts = 0.22;
t=0:0.01:1.5;

Gs = K / (1 + Ts*s); % Transferfunktion PT1-Glied
y = step(Gs,t); 

figure 
plot(data(:,1), data(:,2),'r','LineWidth',2);
hold on;
plot(t,y,'b','LineWidth',2);
grid on;


Autor: Tobias Plüss (hubertus)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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!

Autor: Zitronen Falter (jetztnicht)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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.

Autor: olaf (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Ich staune ueber den mathematischen Level hier. :-)

Kuckt mal hier:

https://de.wikipedia.org/wiki/Polynominterpolation

Habt ihr das nicht in der Schule gelernt?

Für die Faulen unter euch wurde es unter dem namen polyfit in Octave 
eingebaut.

Olaf

Autor: PolynomMaker (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert

Autor: 0815 (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
zuerst ist das Signal zu differenzieren, da es von einem Sprung und 
njchf von ejnem Dirac kommt.
Danach Jnverse fft. Die ANZAHL der Punkte DER inv. fft soll  viel größer 
sein als die Messwerte, (4096 zu 100) da sonst die Messwerte zyklisch 
ausgewertet werden.

Autor: Jan K. (jan_k)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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

: Bearbeitet durch User
Autor: aSma>> (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Servus,

@ Jan K. (jan_k)
ja ich habe gerade dein Ansatz mal soeben ausprobiert, damit kriege ich 
folgende Fkt raus:

G(z)= 4.6z/(z-0.95). Diese passt mit der Übertragungsfkt. von werte.txt 
perfekt überein.

Danke.

Autor: Tobias Plüss (hubertus)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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

Autor: Dirac Impuls (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hallo,

auch von mir ein Danke an Jan!

Trotzdem verstehe ich die Matrix nicht. Das Element unten links 
-y[m+N-1] ergibt für z.B. m=2 -y[N+1]. Das höchste Element in y hat aber 
den Index N.

Autor: Jan K. (jan_k)
Datum:
Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
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.
[y, tt] = stepz(1,[1 -.95 .001 .1],100); % daten generieren
u = ones(size(y)); % eingangssignal
m = 1; % approximation 1. Ordnung
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
p = psi\y(m+1:end); % parameter

figure,
stem(y)
hold all
stepz([0 p(floor(length(p)/2)+1:end)'],[1 p(1:floor(length(p)/2))'])
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.

: Bearbeitet durch User
Autor: Tobias Plüss (hubertus)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
@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.

Autor: aSma>> (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Jan K. schrieb:
> @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.

Servus,
ich schreibe nie etwas im Forum ab. Ist einfach zu Fehleranfällig. Ich 
habe einen Ansatz genommen bei dem man den Exponenten sowohl im Zähler 
als auch im Nenner vorgeben kann. m=0 (Zähler), n=2 (Nenner).

Dies hat auch "Oh Doch (jetztnicht)" vorgeschlagen.

Mit Verlaub, habe ich meiner Formel G(z)= 4.6z/(z-0.95) die Abtastzeit 
vergessen anzugeben: Ta=10ms.

Autor: Jan K. (jan_k)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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? :-)

Autor: Dirac Impuls (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Jan K. schrieb:
> y muss bei
> y(m) starten, nicht bei y(0).

Hallo Jan,

danke nochmal für die Antwort. In dem Fall passt aber das Element oben 
links nicht mehr. Wenn y bei m Anfängt, was ist dann y(m-1)?

Autor: Jan K. (jan_k)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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.

Autor: aSma>> (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Jan K. schrieb:
> Da ist die Namensgebung sehr blöd von mir und beruht auf Faulheit :-/
Deswegen macht man einfach ein Foto. Dieses rumgetippe. Anscheinend 
hasst du nichts bessres zu tun. Such dir mal ne Freundin.

Jan K. schrieb:
> Woher hast du deine Quelle? :-)
Quid pro quo (kirre ausgesprochen). Gibt erstmal eine fehlerfreie 
Version von dir ab.

Dann will ich deine mathematik Künste mal testen:
>>p = psi\y.

Wie berechnet man dies (ohne Matlab)?

Autor: Jan K. (jan_k)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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.

: Bearbeitet durch User
Autor: aSma>> (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Servus,
hier das gewünschte Skript ohne Kommentar:

https://www.tu-ilmenau.de/fileadmin/media/systeman...

Autor: Dirac Impuls (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hallo Jan,

vielen Dank für die Hilfe! Eine Frage hätte ich noch: Wird in dem Buch / 
in den Büchern auch die Systemidentifikation für Mehrgrößensysteme 
beschrieben?

Antwort schreiben

Die Angabe einer E-Mail-Adresse ist freiwillig. Wenn Sie automatisch per E-Mail über Antworten auf Ihren Beitrag informiert werden möchten, melden Sie sich bitte an.

Wichtige Regeln - erst lesen, dann posten!

  • Groß- und Kleinschreibung verwenden
  • Längeren Sourcecode nicht im Text einfügen, sondern als Dateianhang

Formatierung (mehr Informationen...)

  • [c]C-Code[/c]
  • [avrasm]AVR-Assembler-Code[/avrasm]
  • [code]Code in anderen Sprachen, ASCII-Zeichnungen[/code]
  • [math]Formel in LaTeX-Syntax[/math]
  • [[Titel]] - Link zu Artikel
  • Verweis auf anderen Beitrag einfügen: Rechtsklick auf Beitragstitel,
    "Adresse kopieren", und in den Text einfügen




Bild automatisch verkleinern, falls nötig
Bitte das JPG-Format nur für Fotos und Scans verwenden!
Zeichnungen und Screenshots im PNG- oder
GIF-Format hochladen. Siehe Bildformate.
Hinweis: der ursprüngliche Beitrag ist mehr als 6 Monate alt.
Bitte hier nur auf die ursprüngliche Frage antworten,
für neue Fragen einen neuen Beitrag erstellen.

Mit dem Abschicken bestätigst du, die Nutzungsbedingungen anzuerkennen.