| Beschreibung | Der Datensatz X wird spektral analysiert und die spektrale Leistungsdichte (engl.: power spectral density, PSD, das ist das mittlere Amplituden-Quadrat) in Spc gespeichert. xsiz muß mindestens 2*specsiz betragen und specsiz eine ganzzahlige Potenz von 2 sein. Intern wird X in xsiz*specsiz / 2 Segmente unterteilt und über die für die einzelnen Segmente berechneten Spektren gemittelt. Jedes Segment der Länge 2*specsiz ergibt die PSD für specsiz+1 Frequenzen (siehe VF_FFT). Um specsiz als ganzzahlige Potenz von 2 definieren zu können, werden jedoch nur specsiz Punkte in Spc gespeichert. Der letzte Wert, d.h. die Leistungsdichte für die Nyquist-Frequenz fc, wird als Rückgabewert der Funktion behandelt und kann entweder vernachlässigt (durch Aufruf der Funktion wie eine void-Funktion bzw. eine procedure) oder als letztes Element in Spc gespeichert werden durch Aufruf in der Form
Spc[specsiz] = VF_spectrum( Spc, specsiz, X, xsiz, Win );
In diesem Fall muß Spc eine Länge von specsiz+1 haben.
Win ist ein Fenster, das auf die Datensegmente angewendet wird. Seine Größe ist immer 2*specsiz. Innerhalb der VectorLib-Bibliotheken sind drei Funktionen vorhanden, die passende Fenster erzeugen: VF_Welch, VF_Parzen und VF_Hann. Ein Rechteck-Fenster wird über VF_equ1 erhalten. Von seinem Gebrauch wird hier aber abgeraten.
Die Qualität des erhaltenen Spektrums kann auf einfache Weise über das Parseval'sche-Theorem getestet werden (wenn Sie durch Aufruf in der eben skizzierten Form auch den Wert an der Nyquist-Frequenz mit gespeichert haben): 1.0/xsize * VF_ssq( X, xsize ) muß ungefähr gleich VF_sum( Spc, specsiz+1 ) sein. Eine große Abweichung zwischen beiden Ausdrücken deutet auf ein zu großes Abtast-Intervall in X.
Bezüglich spezieller Versionen mit den Präfixen VFp_, VFs_ und
VFl_ vergleiche man Kap. 4.8.. |
| Beispiel C/C++ (für DOS): | #include <VFstd.h>
#include <VFmath.h>
#include <Vgraph.h>
#include <math.h>
#include <conio.h>
void main( void )
{
fVector X, Spc, Win, Time, Freq;
float deltat, fOscill, fNyquist;
ui xsize=4096, specsiz=256;
X = VF_vector( xsize );
Spc = VF_vector( specsiz+1 );
Win = VF_vector( 2*specsiz );
Time = VF_vector( xsize );
Freq = VF_vector( specsiz+1);
deltat = 0.001; /* Abtast-Intervall 1 Millisekunde */
fNyquist = 0.5 / deltat; /* Nyquist-Frequenz = 500 Hz */
fOscill = 100; /* es sei eine 100 Hz-Schwingung */
VF_ramp( Time, xsize, 0, deltat );
VF_ramp( Freq, specsiz+1, 0, fNyquist / specsiz );
VFx_sin( X, Time, xsize, 2*M_PI*fOscill, 0, 1 );
/* Sinus-Schwingung (omega*t) */
VF_cmpC( X, X, xsize, 0.7 );
/* in asymmetrische Rechteckschwingung umwandeln */
VF_Welch( Win, 2*specsiz ); /* oder eine andere Fensterfunktion */
Spc[specsiz] = VF_spectrum( Spc, specsiz, X, xsize, Win );
V_initGraph( "\\BorlandC\\BGI\\" );
/* korrekten BGI-Pfad angeben! */
VF_xyAutoPlot( Freq, Spc, specsiz+1,
PS_SOLID | SY_CROSS, GREEN );
getch(); /* beliebige Taste drücken */
closegraph();
V_nfree( 5, Freq, Time, Win, Spc, X );
} |
| Beispiel Pascal/Delphi: |
uses VFstd, VFmath, Vgraph,...;
const
xsize=4096; specsiz=256;
var
X, Spc, Win, Time, Freq: fVector;
deltat, fOscill, fNyquist: Single;
begin
X := VF_vector( xsize );
Spc := VF_vector( specsiz+1 );
Win := VF_vector( 2*specsiz );
Time := VF_vector( xsize );
Freq := VF_vector( specsiz+1);
deltat := 0.001; (* Abtast-Intervall 1 Millisekunde *)
fNyquist := 0.5 / deltat; (* Nyquist-Frequenz = 500 Hz *)
fOscill := 100; (* es sei eine 100 Hz-Schwingung *)
VF_ramp( Time, xsize, 0, deltat );
VF_ramp( Freq, specsiz+1, 0, fNyquist / specsiz );
VFx_sin( X, Time, xsize, 2*PI*fOscill, 0, 1 );
(* Sinus-Schwingung (omega*t) *)
VF_cmpC( X, X, xsize, 0.7 );
(* in asymmetrische Rechteckschwingung umwandeln *)
VF_Welch( Win, 2*specsiz ); (* oder eine andere Fensterfunktion *)
VF_Pelement( Spc,specsiz )^ :=
VF_spectrum( Spc, specsiz, X, xsize, Win );
VF_xyAutoPlot( Freq, Spc, specsiz+1,
PS_SOLID + SY_CROSS, GREEN );
V_free( Freq ); V_free( Time ); V_free( Win );
V_free( Spc ); V_free( X );
end;
(* Es wird angenommen, daß die Graphik-Routinen bereits initialisiert worden sind (DOS) oder daß das obige Beispiel innerhalb einer Paint- (Pascal für Windows) oder ShowView-Routine (Delphi) steht. *) |
|