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 !