Posted By: pivson (Pijte pivo, je zdrave !) on 'CZprogram' Title: Re: fft v pascalu Date: Wed May 10 17:53:11 2000 Zapomel sem :) {========================================================================== testfft.pas - Don Cross <dcross@intersrv.com> This program is a test/demo for the file 'fourier.pas'. Get the latest version of 'fourier.pas' and 'testfft.pas' at the following URL. http://www.intersrv.com/~dcross/fft.html#pascal NOTE: You may need to modify the const string 'PathToBGI' to point to the correct drive and subdirectory for the BGI drivers on your computer, in order for the graphics to work. --------------- What this program does ------------------------- First, it generates a time signal consisting of a large 200 Hz sine wave added to a small 2000 Hz cosine wave, which is graphed on the screen. (Press ENTER after you are done viewing each graph.) Next, it performs the FFT and graphs the resulting complex frequency samples. Then, it filters out all frequency components above 1000 Hz in the transformed data. Finally, it performs the inverse transform to get a filtered time signal back, and graphs the result. ------------------------ Revision history ------------------------ 1996 December 12 [Don Cross] Added code to test the new procedure Fourier.CalcFrequency. Cleaned up some comments. Added code to preserve the original text mode. 1996 November 17 [Don Cross] Wrote and debugged first version. ==========================================================================} program TestFFT; {$N+,E+} uses Fourier, Graph, Crt; const PathToBGI = '.'; { Change as necessary } function f ( t: double ): double; begin f := sin ( 200 * 2*PI * t ) + 0.2 * cos ( 2000 * 2*PI * t ); end; const NumSamples = 512; { buffer size must be power of 2 } SamplingRate = 22050; { sampling rate in Hz } type Buffer = array [0 .. NumSamples-1] of double; var RealIn, ImagIn, RealOut, ImagOut: Buffer; OutputListingFile: text; i: integer; temp, t, dt: double; procedure Test_CalcFrequency; var yr, yi: double; i: integer; begin { Fill input buffers with random data } for i := 0 to NumSamples-1 do begin RealIn[i] := Random(10000); ImagIn[i] := Random(10000); end; writeln ( OutputListingFile ); writeln ( OutputListingFile, '*** Testing procedure CalcFrequency ***' ); writeln ( OutputListingFile ); fft ( NumSamples, RealIn, ImagIn, RealOut, ImagOut ); for i := 0 to NumSamples-1 do begin CalcFrequency ( NumSamples, i, RealIn, ImagIn, yr, yi ); writeln ( OutputListingFile, i:4, RealOut[i]:15:6, yr:15:6, ImagOut[i]:20:6, yi:15:6 ); end; end; procedure ListData ( var RealData, ImagData: Buffer; comment: string ); var i, yr, yi, prev_yr, prev_yi: word; trash: char; maxAbsValue: double; begin writeln ( OutputListingFile, '*** ', comment, ' ***' ); writeln ( OutputListingFile ); writeln ( OutputListingFile, 'index':20, 'real':20, 'imag':20 ); for i := 1 to NumSamples do begin writeln ( OutputListingFile, i:20, RealData[i]:20:5, ImagData[i]:20:5 ); end; writeln ( OutputListingFile ); writeln ( OutputListingFile, '------------------------------------------------------------------------' ); writeln ( OutputListingFile ); maxAbsValue := 0.0; for i := 0 to NumSamples-1 do begin if abs(RealData[i]) > maxAbsValue then maxAbsValue := abs(RealData[i]); if abs(ImagData[i]) > maxAbsValue then maxAbsValue := abs(ImagData[i]); end; for i := 0 to NumSamples-1 do begin yr := Trunc ( GetMaxY * (1 - (RealData[i] / maxAbsValue + 1)/2) ); yi := Trunc ( GetMaxY * (1 - (ImagData[i] / maxAbsValue + 1)/2) ); if i > 0 then begin SetColor ( LIGHTRED ); Line ( i-1, prev_yr, i, yr ); SetColor ( LIGHTGREEN ); Line ( i-1, prev_yi, i, yi ); end; prev_yr := yr; prev_yi := yi; end; sound (800); delay (100); nosound; trash := ReadKey; (* Pause *) ClearDevice; end; var GraphDriver, GraphMode, StartupTextMode: integer; FreqIndex: word; begin StartupTextMode := LastMode; assign ( OutputListingFile, 'fftout.txt' ); rewrite ( OutputListingFile ); GraphDriver := VGA; GraphMode := VGAHI; InitGraph ( GraphDriver, GraphMode, PathToBGI ); dt := 1.0 / SamplingRate; t := 0.0; for i := 0 to NumSamples-1 do begin RealIn[i] := f(t); ImagIn[i] := 0.0; t := t + dt; end; ListData ( RealIn, ImagIn, 'Time domain data before transform' ); fft ( NumSamples, RealIn, ImagIn, RealOut, ImagOut ); ListData ( RealOut, ImagOut, 'Frequency domain data after transform' ); (* Filter out everything above 1000 Hz (low-pass) *) FreqIndex := Trunc ( 1000.0 * NumSamples / SamplingRate ); for i := 0 to NumSamples-1 do begin if ((i > FreqIndex) and (i < NumSamples DIV 2)) or ((i >= NumSamples DIV 2) and (i < NumSamples-FreqIndex)) then begin RealOut[i] := 0.0; ImagOut[i] := 0.0; end; end; ifft ( NumSamples, RealOut, ImagOut, RealIn, ImagIn ); ListData ( RealIn, ImagIn, 'Time domain data after inverse transform' ); Test_CalcFrequency; close ( OutputListingFile ); CloseGraph; TextMode ( StartupTextMode ); end. {--- end of file testfft.pas ---} Pivson I a posledni, z bozi vule pivar http://pulse.mute.cz http://it3.mute.cz A co budou delat cesi ??? Deme na pivo !