#!/usr/bin/octave -q

# Compare the results of kiss_fft with those calculated by octave's fft routines

nfft=1024;
showdiff=0;

for test_case =1:3
    printf( '====================================\n');
    if test_case == 1,
        prec = 'short'
        prog = './kiss_fft_s'
        scale = nfft
    elseif test_case  == 2,
        prec = 'float'
        prog = './kiss_fft_f'
        scale = 1
    else
        prec = 'double'
        prog = './kiss_fft_d'
        scale = 1
    endif

    # create the input
    input = rand(1,nfft) .* exp(2*pi*j*rand(1,nfft));
    input = floor( real(input)*32767 ) + j*floor(imag(input)*32767 );

    # write the input
    f = fopen( "fftin.dat", "w", "native");
    to_write = [real(input);imag(input) ];
    fwrite(f,to_write,prec);
    fclose(f);

    cmd = sprintf('%s %d < fftin.dat > fftout.dat',prog,nfft)
    system(cmd);

    # read the output
    f = fopen("fftout.dat","r", "native");
    from_read = fread(f,Inf,prec)';
    output = (from_read(1:2:nfft*2) + j*from_read(2:2:nfft*2))        ;

    comp = fft(input)/scale;
    if test_case == 1,
        comp = floor( real(comp) +.5) + j*floor( imag(comp) +.5);
    endif

    if showdiff,
        x=linspace(0,2*pi,nfft);
        figure(1);
        plot( x, real(output), x, real(comp) );
        figure(2);
        plot( x, imag(output), x, imag(comp) );
        figure(3);
        plot( x, abs(output-comp) );
        pause
    else
        diff = output - comp;
        noise_pow =  sum( conj(diff).*diff );
        sig_pow = sum( conj(comp).*comp );
        snr = 10*log10( sig_pow / noise_pow )

        avg_scale = mean(  abs(output) ./ abs(comp) )
    endif
endfor