Impulsantwort aus Frequenzgang durch IFFT in Matlab

[Janek U. (patalmqx)]

Guten Abend liebe Götter des DSP!

Ich sitze wieder einmal vor Matlab und habe diesmal einen Frequenzgang 
eines Balkenschwingers (SISO, Eingang: Kraftaufnehmer am 
elektrodynamischen Shaker, Ausgang: 1 Beschleunigungssensor auf dem 
Balken) mittels der Welch Methode ermittelt.
Dieser stimmt auch mit dem Frequenzgang einer kommerziellen 
Schwingungsanalyse Software überein, soweit ich das sehen kann.

Der Frequenzgang ist in der frf.png zu sehen.
Der dazugehörige Frequenzvektor ist: fv=0:0.25:1599.75 [Hz] (Länge 
6400). Die DC Komponente der frf ist reell.

Nun kann man sich die Impulsantwort des Systems ja meines Wissens nach 
mit einer ifft berechnen. Dazu habe ich das Spektrum folgendermaßen 
erweitert (einseitig zu zweiseitig):

h1n=[0.5*h1 h1(1) 0.5*fliplr(conj(h1(2:end)))];

Der Frequenzvektor wird ebenfalls 'verlängert': (Länge nun 12800)

fvp=0:df:2*fv(end)+df;
fvp=0:0.25:3199.75;

Ifft dieser FRF h1n ergibt, wie erwartet, ein reelles Signal. Am Anfang 
besitzt es jedoch mit fast -9 einen sehr starken Ausschlag, der sofort 
abnimmt. Den Anstieg gegen Ende der IRF deute ich mal als mathematische 
Ungenauigkeiten/Messungenauigkeiten/Nichtlinearitäten..?

Habe ich etwas übersehen, oder ist diese Impulsantwort für dieses Modell 
richtig? Oder ist der starke Ausschlag eine Folge der Annahme als 
LZI-System?  Bei Bedarf kann ich auch den Code plus Daten zur Verfügung 
stellen!

Für Tips wäre ich sehr dankbar, da Beispiele zu diesem Problem sehr rar 
gesät sind!

Danke, patalmqx

[Detlef _. (detlef_a)]

Bei 1600Hz hast Du noch fette Frequenzkomponenten. Du solltest für hohe 
Frequenzen keine Spektralanteile mehr haben. Im Zeitbereich heißt das, 
dass Nyquist erfüllt sein muß.

Der 'mittlere' Wert, bei Dir h1(1) ist der Nyquistanteil, den setzt Du 
0.

>>Den Anstieg gegen Ende der IRF deute ich mal als mathematische
Ungenauigkeiten/Messungenauigkeiten/Nichtlinearitäten..?

Glaub ich nicht, das liegt an was anderem, vllt. wie die Daten 
gefenstert/gemittelt werden. Was ist mit den Phasen, ich dachte Welch 
liefert ne PSD? Wenn Du alle Phasen Null setzt kommt genau der fette 
peak bei t=0 raus.

Cheers
Detlef

[Janek U. (patalmqx)]

Danke erstmal für die Antwort!
Zu den Spektralanteilen: Ich habe die Sensor-Position auf dem Balken 
variiert und damit mehrere FRFs aufgenommen. Wenn nun bei 1550 Hertz 
eine Eigenfrequenz des Balkens sitzt, muss ich die doch auch in der 
Amplitude sehen können, oder nicht?

Mit Nyquist meinst Du wahrscheinlich das Abtasttheorem? Ich denke, dass 
ich das nicht verletzt habe. Ich habe in meiner Software die maximale 
Frequenz im Spektrum auf 1,6kHz eingestellt, worauf es einen AA Filter 
einsetzt. Die Rohdaten (Kraft, Beschleunigung) zeigen auch, dass mit 
4,096 kHz abgetastet wurde, also dem 2.56 fachen von 1,6 kHz..für mich 
ein Signal, dass ein AA Filter eingesetzt wird. Bei der hochpreisigen 
Software würde mich ein fehlender AA Filter auch wundern. Ebenfalls ist 
ein 0,7 Hz Hochpass eingestellt.

Zum Nyquistanteil: Wieso setze ich den gleich 0? Ich hatte mit meiner 
"Methode" mal ein einfaches Beispiel (f(t)->F(f) einseitig->F(f) 
zweiseitig -> (IFFT) f(t)) programmiert. Dieses hatte funktioniert. Hier 
der Code:

clear;close all;
%23.10.13
%FRF2IRF test

n=1024; %Länge der Daten/FFT
fs=30;  %Abtastfrequenz
t=n/fs; %Gesamtzeit
dt=t/n; %Zeitinkrement
df=1/(t-dt);  %Frequenzinkrement im Spektrum

tv=0:dt:t-dt; %Zeitvektor
fv=0:df:fs; %Frequenzvektor

freq1=1;%Hz
freq2=0.5;%Hz
for j=1:length(tv)
    tj=j*dt; %aktuelle Zeit
    u(j)=0.5*sin(2*pi*freq1*tj)+cos(2*pi*freq2*tj); %Zeitsignal
end
plot(tv,u); hold on;

%Spektrum erzeugen mit FFT
sp=(1/length(u))*fft(u);

%Originalsignal durch IFFT
u_selbst=length(u)*ifft(sp);
%plot(tv,u_selbst,'r--');xlabel('t[s]');ylabel('u(t)');

%Einseitiges Spektrumn
lh=length(u)/2; %halbe Länge
hsp=2*sp(1:lh); %doppelte Amplitude

%Einseitiges Spektrum zu zweiseitigem Spektrum zusammensetzen
sp_s=0.5*[hsp hsp(1) fliplr(conj(hsp(2:end)))];

figure();
plot(fv,abs(sp),fv,abs(sp_s),'r--');title('zweiseitiges Spektrum aus FFT und aus einseitigem Spektrum zusammengesetzt (gestrichelt)');

fu_s=length(u)*ifft(sp_s); %Zeitsignal durch IFFT des zweiseitgen Spektrums wieder herstellen
figure();
plot(tv,fu_s,tv,u_selbst,'k')
hold on;
plot(tv,u,'r--');title('Originales Signal (rot) und aus IFFT erzeugtes Signal (schwarz)');

Angehangen sind die FRF von einer anderen Position am Balken aufgenommen 
plus PHase. Die rote Kurve ist die IRF aus meinem FFT Programm (noch 
ohne Frequenzachse) und die grüne IRF aus dem Programm. Die Kurven sehen 
sich doch schon irgendwie ähnlich..
Wieso ist die Einheit der IRF des Programms eigentlich 
(m/s^2)/N/s=1/(kg*s)? Eigenartig, ich hätte erstmal m/s^2 vermutet.

In meinem Messungen gab es bei 1-1,2 kHz meist Oszillationen in den FRFs 
sowie eine schlechte Kohärenz, wahrscheinlich sind die Daten von 
schlechter Qualität. Das könnte zumindest die Abweichung zur 
Programm-IRF erklären, bei dieser FRF Messung traten nämlich keine 
Oszillationen auf.

Ich muss mich damit noch eingehender beschäftigen, besser spamme ich 
nicht das Forum zu ;)

[patalmqx (Gast)]

Für die Nachwelt: Das oben gezeigte Phänomen kann bei nichtlinearen 
Systemen auftreten. Neben einem linearen System und möglichst guter 
Datenqualität kann man das je nachdem noch mit der 
Hilbert-Transformation behandeln (Anders Brandt, "Noise and Vibration 
Analysis" behandelt das z.B.).

Die Lösung scheint aber vor allem zu sein, dass man für die irf nur die 
Hälfte der Werte aus der IFFT des doppelseitigen Spektrums nehmen darf, 
um die IRF zu erhalten (wenn ich mich nicht irre). So komme ich mit dem 
Zeitvektor wieder genau auf meine Messdauer T=4s und erhalte eine IRF 
die auch zu einem kausalen System passt.

So, gute Tat getan.