FFT.c File Reference

Fast Fourier Trasform Functions. More...

#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"

Include dependency graph for FFT.c:

Go to the source code of this file.

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

Functions
int	CheckSpikes (int mod, double *tmp, int npix)
int	FilterData (unsigned short *mat, int horpix, int verpix, int filw, int bp)
void	four1 (double data[], unsigned long nn, int isign)
void	realft (double data[], unsigned long n, int isign)
int	RemoveSpikes (int mod, unsigned short *mat, int horpix, int verpix, int blindpix)
	Removes Spikes on a measure. More...
int	SmoothingF (double *tmp, int npix, int pts)
	Smoothing of 1D array. More...

Variables
doas	DOAS

Detailed Description

Fast Fourier Trasform Functions.

Author

Daniele Bortoli

Definition in file FFT.c.

Macro Definition Documentation

SWAP

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

Definition at line 24 of file FFT.c.

Referenced by four1().

Function Documentation

CheckSpikes()

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

Definition at line 348 of file FFT.c.

Referenced by RemoveSpikes().

{

    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;
}

Here is the caller graph for this function:

FilterData()

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

Definition at line 250 of file FFT.c.

References IDX, and SmoothingF().

Referenced by BuildLogRatio(), DB_Column(), and DrawAreaPlotCB().

{


    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;

}

Here is the call graph for this function:

Here is the caller graph for this function:

four1()

void four1	(	double	data[],
		unsigned long	nn,
		int	isign
	)

Definition at line 26 of file FFT.c.

References SWAP.

Referenced by realft().

{
    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;
    }
}

Here is the caller graph for this function:

realft()

void realft	(	double	data[],
		unsigned long	n,
		int	isign
	)

Definition at line 95 of file FFT.c.

References c2, and four1().

Referenced by SmoothingF().

{
    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.

    }
}

Here is the call graph for this function:

Here is the caller graph for this function:

RemoveSpikes()

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

Removes Spikes on a measure.

Parameters

mod	NOT USED (=0)
*mat	source buffer
horpix	horizontal pixels
verpix	vertical pixels
blindpix	blind pixels

Definition at line 369 of file FFT.c.

References CheckSpikes(), and IDX.

Referenced by DrawAreaPlotCB().

{


    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;

}

Here is the call graph for this function:

Here is the caller graph for this function:

SmoothingF()

int SmoothingF	(	double *	tmp,
		int	npix,
		int	pts
	)

Smoothing of 1D array.

Parameters

*tmp	source buffer
npix	number of elements of the array
pts	smoothing window

Definition at line 151 of file FFT.c.

References MMessageDialog(), and realft().

Referenced by FilterData(), SmoothData(), and SmoothFloatData().

{
    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;
}

Here is the call graph for this function:

Here is the caller graph for this function:

Variable Documentation

DOAS

doas DOAS

Definition at line 15 of file FFT.c.