#!/usr/bin/octave -q

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

nfft=1024
showdiff=0;
use_inverse_fft = 1

infile = 'fftin.dat'
outfile = 'fftout.dat'
for test_case =1:3
    printf( '====================================\n');
    if test_case == 1,
        prec = 'short'
        prog = './kiss_fft_s'
    elseif test_case  == 2,
        prec = 'float'
        prog = './kiss_fft_f'
    else
        prec = 'double'
        prog = './kiss_fft_d'
    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( infile, "w", "native");
    to_write = [real(input);imag(input) ];
    fwrite(f,to_write,prec);
    fclose(f);

    if use_inverse_fft,
        cmd = sprintf('%s -s -i %d < %s > %s',prog,nfft,infile,outfile)
        comp = ifft(input);
    else        
        cmd = sprintf('%s -s %d < %s > %s',prog,nfft,infile,outfile)
        comp = fft(input);
    endif
    system(cmd);

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

    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) )
    var_scale = var(  abs(output) ./ abs(comp) )
    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
    endif
endfor