Grundschwingung in MATLAB berechnen

[Alex (Gast)]

Hi,

ich habe den Spannungsverlauf an einer Asynchronmaschine mit dSpace 
aufgenommen und möchte jetzt gerne die Grundschwingung / 1. Harmonische 
in MATLAB berechnen, da ich für weitere Arbeiten den 
Grundschwingungseffektivwert usw. berechnen will.


%Speichern der Messwerte in Vektor
t=mw(:,1);
u1=mw(:,2);

Die beiden Vektoren habe eine Länge von 1024000, wobei genau 20 Perioden 
gespeichert wurden (2 s Aufnahmezeit bei 10 Hz).
Wie kann man nun die 1. Harmonische ermitteln? Ist das richtig, dass ich 
zuerst einmal den Datenvektor auf 1 Periode beschneiden sollte?

% Daten auf 1 Periode zuschneiden
u1 = u1(1:length(u1)/20);
t = t(1:length(t)/20);

Ist meine Berechnung der Grundschwingung korrekt?

% Berechnung der Grundschwingung
N = length(u1);
 for n=1:N
   u1gs(n) = u1(n)*exp(-1i*((2*pi)/N)*n);
 end


Ich würde euch bitten mir da zu helfen.

[Yalu X. (yalu) (Moderator)]

Alex schrieb:
> Ist das richtig, dass ich zuerst einmal den Datenvektor auf 1 Periode
> beschneiden sollte?

Das kannst du tun, must es aber nicht. Wichtig ist nur, dass du eine
ganzzahlige Anzahl von Perioden bearbeitest. Mehr Perioden bedeutet mehr
Rechenzeit, aber auch eine bessere Unterdrückung des Messrauschens.

> u1gs(n) = u1(n)*exp(-1i*((2*pi)/N)*n);

Stimmt, allerdings ist das noch nicht die Grundschwingung selbst. Von
u1gs musst du noch den Mittelwert bilden, diesen mit 2 multiplizieren,
und schon hast du die komplexe Amplitude der Grundschwingung. Diesen
kannst du dann in Betrag und Phase zerlegen. Für die Mittelwertbildung
gibt es in Matlab sicher eine fertige Funktion, so dass du das Ganze
ohne (ineffiziente) Schleife hinschreiben kannst.

[Alex (Gast)]

Ich habe jetzt die Zeilen so geändert (noch ineffizient, aber für's 
erste für mich nicht wichtig), wie du es mir vorgeschlagen hast:

N = length(u1);
 for n=1:N
   u1gs(n) = 2/N*(u1(n)*exp(-1i*((2*pi)/N)*n));
 end

Im Anhang siehst du die Messdaten in einem Plot, die ich auswerte. Ich 
würde jetzt gerne die Grundschwingung als Zeitfunktion in den Plot 
einzeichnen. Ist das mit der komplexen Zahl, die ich bei u1gs 
rausbekomme, überhaupt möglich?
Wenn ich Betrag und Phase mit "abs" und "angle" ausrechne, wie komme ich 
dann auf die Zeitfunktion für die Grundschwingung?

[jetzt (Gast)]

Wenn die Frequenz bekannt ist kann man gleich eine Vereinfachte fourier 
mit nur einer Frequenz machen. Die Grundschwingung ist dann :

realteil = 1/T integral f(t) sin ft dt
Imaginaer  = 1/T integral f(t) cos ft dt

Fuer die Oberwellen N = 2..M dann :

realteil = 1/T integral f(t) sin Nft dt
Imaginaer  = 1/T integral f(t) cos Nft dt

[Yalu X. (yalu) (Moderator)]

Alex schrieb:
> Wenn ich Betrag und Phase mit "abs" und "angle" ausrechne, wie komme ich
> dann auf die Zeitfunktion für die Grundschwingung?

u(t) = abs(a) * cos(omega * t + angle(a))

a ist dabei die zuvor berechnete komplexe Amplitude, omega die
(bekannte) Kreisfrequenz der Grundschwingung.

[Alex (Gast)]

Ich habe jetzt die Grundschwingung mit der Formel von dir berechnet.

u1 = u1(1:length(u1)/20);
t = t(1:length(t)/20);

N = length(u1);
for n=1:N
   u1gs(n) = 2/N*(u1(n)*exp(-1i*((2*pi)/N)*n));
   i1gs(n) = 2/N*(i1(n)*exp(-1i*((2*pi)/N)*n));
end

for n=1:length(abs(u1gs))
   u1gs_t(n) = abs(u1gs(n))*cos(2*pi*10*t(n)+angle(u1gs(n)));
   i1gs_t(n) = abs(i1gs(n))*cos(2*pi*10*t(n)+angle(i1gs(n)));
end

Interessanterweise habe ich aber nur so eine geringe Amplitude im 
Millibereich erhalten. Wenn ich mir aber die FFT von u1 anschaue, so 
habe ich aber einen viel, viel größeren Peak bei dieser Frequenz.
Da muss irgendwas falsch sein bei der Berechnung der Grundschwingung.

[Yalu X. (yalu) (Moderator)]

Alex schrieb:
> Da muss irgendwas falsch sein bei der Berechnung der Grundschwingung.

Ja, da fehlt die Mittelwertbildung.

Weil ich Matlab weder habe noch kann, habe ich die Berechnung mal in
Python programmiert, was du sicher ganz leicht in Matlab übersetzen
kannst:

from math import sin, cos, pi, radians, degrees
from cmath import exp, phase
from sys import stderr

f     = 10              # Frequenz
T     = 1 / f           # Periodendauer
omega = 2 * pi * f      # Kreisfrequenz
n     = 1024000 // 20   # Anzahl Samples in 1 Periode
dt    = T / n           # Zeit zwischen zwei aufeinanderfolgenden Samples

# Funktion für die Messreihe mit 7 Harmonischen
# Die Grundschwingung hat den Betrag 3.0 und die Phasenverschiebung 72°

def f(t):
  return   3.0 * cos(1 * omega * t + radians(72)) \
         - 2.4 * cos(2 * omega * t + radians(10)) \
         + 1.9 * cos(3 * omega * t - radians(20)) \
         + 1.7 * cos(4 * omega * t + radians(30)) \
         - 1.5 * cos(5 * omega * t - radians(40)) \
         - 1.3 * cos(6 * omega * t + radians(50)) \
         + 1.1 * cos(7 * omega * t + radians(60))

# gen_t1_u1 generiert eine Beispielsmessreihe aus der obigen Funktion f

def gen_t1_u1():
  t1  = []
  u1 = []
  for i in range(n):   # 1 Periode
    t = i * dt
    t1.append(t)
    u1.append(f(t))
  return t1, u1

t1, u1 = gen_t1_u1()

# Berechnung der Grundschwingung

n = len(u1)
summe = 0
for i in range(n):
  summe += u1[i] * exp(-1j * (2*pi*i/n))
mittelwert = summe / n
amplitude = 2 * mittelwert

betrag = abs(amplitude)
phi = phase(amplitude)

print('betrag = %.2f  phi = %.2f°' % (betrag, degrees(phi)))

# Ausgabe für Plot

plot = open('data', 'w')
for i in range(n):
  t = t1[i]
  gs = betrag * cos(omega * t + phi)
  plot.write('%f %f %f\n' % (t, u1[i], gs))
plot.close()

Die Ausgabe des Programms:

betrag = 3.00  phi = 72.00°

Diese Werte stimmen genau mit den Parametern des ersten Summanden in der
Funktion f überein. Im Anhang sind u1 und die berechnete Grundschwingung
über die Zeit geplottet.