FFT.c

Go to the documentation of this file.

#include <windows.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <malloc.h>
#include <math.h>
#include <direct.h>
#include "mgui.h"
#include "DAS_Spatram.h"
#include "dcl.h"
#include "dil.h"
#include "bil.h"
#include "DOAS.h"

doas DOAS;

#define SWAP(a,b) temp=(a);(a)=(b);(b)=temp

void four1(double data[], unsigned long nn, int isign)
/*
Replaces data[1..2*nn] by its discrete Fourier transform, if isign is input as 1; or replaces
data[1..2*nn] by nn times its inverse discrete Fourier transform, if isign is input as - 1.
data is a complex array of length nn or, equivalently, a real array of length 2*nn. nn MUST
be an integer power of 2 (this is not checked for!).
*/
{
    unsigned long n,mmax,m,j,istep,i;
    double wtemp,wr,wpr,wpi,wi,theta; //Double precision for the trigonomet-ricrecurrences
    double temp,tempi;
    n=nn << 1;
    j=1;

    for (i=1;i<=n;i+=2)
    {  //This is the bit-reversal section of the routine. 
        if (j > i)
        {
            SWAP(data[j],data[i]); //Exchange the two complex numbers.
            SWAP(data[j+1],data[i+1]);
        }
        m=nn;
        while (m >= 2 && j > m)
        {
            j -=m;
            m >>= 1;
        }
        j +=m;
    }
//Here begins the Danielson-Lanczos section of the routine.
    mmax=2;
    while (n > mmax) 
    { //Outer loop executed log 2 nn times.
        istep=mmax << 1;
        theta=isign*(6.28318530717959/mmax); //Initialize the trigonometric recurrence.
        wtemp=sin(0.5*theta);
        wpr = -2.0*wtemp*wtemp;
        wpi=sin(theta);
        wr=1.0;
        wi=0.0;
        
        for (m=1;m<=mmax;m+=2)
        { //Here are the two nested inner loops.
            for (i=m;i<=n;i+=istep) 
            {
                j=i+mmax; //This is the Danielson-Lanczos formula:
                temp=wr*data[j]-wi*data[j+1];
                tempi=wr*data[j+1]+wi*data[j];
                data[j]=data[i]-temp;
                data[j+1]=data[i+1]-tempi;
                data[i] += temp;
                data[i+1] += tempi;
            }
            wr=(wtemp=wr)*wpr-wi*wpi+wr; //Trigonometric recurrence.
            wi=wi*wpr+wtemp*wpi+wi;
        }
        mmax=istep;
    }
}

/*
Calculates the Fourier transform of a set of n real-valued data points. Replaces this data (which
is stored in array data[1..n]) by the positive frequency half of its complex Fourier transform.
The real-valued .rst and last components of the complex transform are returned as elements
data[1] and data[2], respectively. n must be a power of 2. This routine also calculates the
inverse transform of a complex data array if it is the transform of real data. (Result in this case
must be multiplied by 2/n.)
*/

void realft(double data[], unsigned long n, int isign)
{
    void four1(double data[], unsigned long nn, int isign);
    unsigned long i,i1,i2,i3,i4,np3;
    double c1=0.5,c2,h1r,h1i,h2r,h2i;
    double wr,wi,wpr,wpi,wtemp,theta; //Double precision for the trigonometric recurrences.
    //theta=3.141592653589793/(double) (n>>1);// Initialize the recurrence.
    theta=3.141592653589793/(double) (n);// Initialize the recurrence.
    if (isign == 1)
    {
        c2 = -0.5;
        //four1(data ,n>>1,1); //The forward transform is here.
        four1(data ,n ,1); //The forward transform is here.
    }
    else
    {
        c2=0.5; //Otherwise set up for an inverse transform.
        theta = -theta;
    }
    wtemp=sin(0.5*theta);
    wpr = -2.0*wtemp*wtemp;
    wpi=sin(theta);
    wr=1.0+wpr;
    wi=wpi;
    //np3=n+3;
    np3=2 *n + 3;
    //for (i=2;i<=(n>>2);i++)           
    for (i=2;i<=(n/2)+1;i++)            
    { //Case i=1 done separately below.
        i4=1+(i3=np3-(i2=1+(i1=i+i-1)));
        h1r=c1*(data[i1]+data[i3]); //The two separate transforms are separated out of data.
        h1i=c1*(data[i2]-data[i4]);
        h2r = -c2*(data[i2]+data[i4]);
        h2i=c2*(data[i1]-data[i3]);
        data[i1]=h1r+wr*h2r-wi*h2i; //Here they are recombined to form the true transform of the original real data.
        data[i2]=h1i+wr*h2i+wi*h2r;
        data[i3]=h1r-wr*h2r+wi*h2i;
        data[i4] = -h1i+wr*h2i+wi*h2r;
        wr=(wtemp=wr)*wpr-wi*wpi+wr; //The recurrence.
        wi=wi*wpr+wtemp*wpi+wi;
    }
    if (isign == 1)
    {
        data[1] = (h1r=data[1])+data[2]; //Squeeze the .rst and last data together to get them all within the original array.
        data[2] = h1r-data[2];

    }
    else
    {
        data[1]=c1*((h1r=data[1])+data[2]);
        data[2]=c1*(h1r-data[2]);
        //four1(data ,n>>1,-1); //This is the inverse transform for the case isign=-1.
        four1(data , n, -1); //This is the inverse transform for the case isign=-1.

    }
}
int SmoothingF(double *tmp, int npix, int pts)
{
    long double a, kost;
    int NMin = 0;
    //int mmax = 4096;
    int mmax = 8192;
    int    M = 2, i, k;
    unsigned long MO2;
    long double y0, yn, FAC;
    long double rn1, dummy, dummy1;
    
    NMin = npix + 2 * pts;
    do
    { 
       M = 2 * M;   // find the next power of 2
    
    }while(M < NMin);
    
    if(M > mmax)
    {
        MMessageDialog("Smoothing Error", "MMax too small", "ok", NULL);
        return 1;   
    } 
        
    kost = pow( (double)pts / (double)M, (double)2.0); //useful constant
    
    
    y0 = tmp[1];
    yn = tmp[npix ];
    rn1 = (double) (1.0 / (double)(npix - 1));
/* 
    //remove the linear trend
    for(i = 0; i< npix;i++)
        tmp[i] = tmp[i] - rn1 * (y0 * (double)(npix - i) + yn * (double)(i - 1));
*/
    //remove the linear trend
    for(i = 1; i<= npix;i++)
        tmp[i] = tmp[i ] - rn1 * (y0 * (double)(npix - i) + yn * (double)(i - 1));
    
    if(npix+1  <= M)
    {
        for(i = npix+1; i<M;i++)
            tmp[i] = 0;
    }

    MO2 = M / 2;
    realft(tmp , MO2, 1); //fourier transform
      
    tmp[1] = tmp[1] / MO2;
    FAC = 1; // square windows 
      
    for (i=1; i <= (int) (MO2 - 1);i++)
    {
        k = 2 * i + 1;     
        
        if ((FAC > 0) | (FAC <0))
            FAC = 0;
       
        a = (1 - kost * pow((double)i, (double)2)) / MO2;
        //a = (1 - kost * pow((0.25 * i),2)) / MO2;
        //a = (1 - pow((0.25 * i),2)) / MO2;
        //a = 1 - pow( (((double)i - MO2/2)/( MO2/2)),2 );
        
        if(a > FAC) 
        {
            FAC = a;
            tmp[k] = FAC * tmp[k];
            tmp[k + 1] = FAC * tmp[k + 1];
        }
        else
        {
            tmp[k] = 0;
            tmp[k + 1] = 0;
        }
    }
    FAC = 0;
    dummy = pow(((double)0.25 * (double)pts),(double)2);
    dummy1= (0.25 * pts) * (0.25 * pts);
    a = (1 - pow(((double)0.25 * (double)pts),(double)2)) / MO2;
    if (a > 0)
        FAC = a;

    tmp[2] = FAC * tmp[2];
       
    realft(tmp , MO2, -1); //inverse fourier transform 
 /*     
    for (i=0;i<npix;i++)
        //Re-Insert the linear trend
        tmp[i] = tmp[i]  + rn1 * (y0 * (double)(npix - i) + yn * (double)(i - 1));
*/
    for (i=1;i<=npix;i++)
        //Re-Insert the linear trend
        tmp[i] = (tmp[i]  + rn1 * (y0 * (double)(npix - i) + yn * (double)(i - 1)));
        
    
    return 0;
}


int FilterData(unsigned short *mat, int horpix, int verpix, int filw, int bp)
{


    int x, y, n ;
    FILE *fd;
    double *avect;
    unsigned long eln = (unsigned long)(2 * (horpix - bp + 1))  ;

    
    fd = fopen("FilterUS_In.tmp", "w");
    if(fd < 0) return 1;
//  for(x = 0; x < horpix + 1; x++)
    for(x = 0; x < horpix ; x++)
    {
        if ((x <= horpix - bp ) & (x > 1))
        {
            fprintf(fd, "%04u ", x);
    
            for(y = 0; y < verpix; y++)
            {   
                fprintf(fd, "%05u ", mat[y * IDX + x]);
            }
            fprintf(fd, "\n", NULL);
        }
    }
    fprintf(fd, "\n", NULL);
    fclose(fd);



    if(filw <= 512)
        avect = (double *) calloc( eln, sizeof(double));
    else if((filw > 512) & (filw <1024))
        avect = (double *) calloc( 2*eln, sizeof(double));
    else if((filw >= 1024) & (filw <2048))
        avect = (double *) calloc( 4*eln, sizeof(double));


    if( avect == NULL )
        return 1;

    for(y = 0; y < verpix; y++)
    {
        n=1;
        for(x = 0; x < horpix; x++)
        {
            if ((x < horpix - bp) & (x > 1))
                {
                    avect[n] = (double)mat[y * (horpix) + x];
                    n++;
                }
            
        }

        SmoothingF(avect, n - 1 , filw);
        
        for(x = 0; x < horpix; x++)
        {
            if ((x < horpix - bp) & (x > 1))
            {
                mat[y * horpix + x] = (unsigned int)avect[x - 1] ;
            }

        }
    }
    
    free(avect);



    fd = fopen("FilterUS_Out.tmp", "w");
    if(fd < 0) return 1;

//  for(x = 0; x < horpix + 1; x++)
    for(x = 0; x < horpix ; x++)
    {
        if ((x <= horpix - bp ) & (x > 1))
        {
            fprintf(fd, "%04u ", x);
    
            for(y = 0; y < verpix; y++)
            {   
                fprintf(fd, "%05u ", mat[y * IDX + x]);
            }
            fprintf(fd, "\n", NULL);
        }
    }
    fprintf(fd, "\n", NULL);
    fclose(fd);

    return 0;

}




int CheckSpikes(int mod, double *tmp, int npix)
{

    double deltapoint = 0;
    double deltapoint1 = 0;

    int i;

    for (i = 1;i<npix - 1;i++)
    {
        deltapoint = fabs(   (double)( (tmp[i+1] - tmp[i])/60000*100)     );
        deltapoint1 = fabs(   (double)( (tmp[i+1] - tmp[i+2])/60000*100)     );
        if( (deltapoint > 1.6) && (deltapoint1 > 1.6)) //1.6 = 1000/60000*100
            tmp[i+1] = tmp[i];

    }

    return 0;
}


int RemoveSpikes(int mod, unsigned short *mat, int horpix, int verpix, int blindpix)
{


    int x, y, n ;
    double *avect;
    unsigned long eln = (unsigned long)(2 * (horpix - blindpix + 1))  ;
    FILE *fd;
    
    fd = fopen("Spikes.tmp", "w");
    if(fd < 0) return 1;
//  for(x = 0; x < horpix + 1; x++)
    for(x = 0; x < horpix ; x++)
    {
        if ((x <= horpix - blindpix ) & (x > 1))
        {
            fprintf(fd, "%04u ", x);
    
            for(y = 0; y < verpix; y++)
            {   
                fprintf(fd, "%05u ", mat[y * IDX + x]);
            }
            fprintf(fd, "\n", NULL);
        }
    }
    fprintf(fd, "\n", NULL);
    fclose(fd);



    avect = (double *) calloc( eln, sizeof(double));


    if( avect == NULL )
        return 1;

    for(y = 0; y < verpix; y++)
    {
        n=1;
        for(x = 0; x < horpix; x++)
        {
            if ((x < horpix - blindpix) & (x > 1))
                {
                    avect[n] = (double)mat[y * (horpix) + x];
                    n++;
                }
            
        }

        CheckSpikes(0, avect, n - 1);
        
        for(x = 0; x < horpix; x++)
        {
            if ((x < horpix - blindpix) & (x > 1))
            {
                mat[y * horpix + x] = (unsigned int)avect[x - 1] ;
            }

        }
    }
    
    free(avect);



    fd = fopen("NoSpikes.tmp", "w");
    if(fd < 0) return 1;

//  for(x = 0; x < horpix + 1; x++)
    for(x = 0; x < horpix ; x++)
    {
        if ((x <= horpix - blindpix ) & (x > 1))
        {
            fprintf(fd, "%04u ", x);
    
            for(y = 0; y < verpix; y++)
            {   
                fprintf(fd, "%05u ", mat[y * IDX + x]);
            }
            fprintf(fd, "\n", NULL);
        }
    }
    fprintf(fd, "\n", NULL);
    fclose(fd);

    return 0;

}