Pełna transformata Fouriera z filtrem dolnoprzepustowym – C/C++

//*****PELNA TRANSFORMATA FOURIERA Z FILTREM DOLNOPRZEPUSTOWYM*******************************
// ©2015 Patryk Grądys
#include 
#include 
#include 

const int len = 1024;

struct CMPL //struktura liczby zespolonej
{
 double Re;
 double Im;
};

struct Fourier
{
 CMPL DFT[len];
 CMPL iDFT[len];
}; 

int main()
{
 int i, k, n, syg[len];
 double x;

 struct Fourier *p;  
 p  = (Fourier*)malloc(sizeof(Fourier));
 
 FILE *pFile;
 
 for(n = 0; n < len; n++)
 {
  syg[n] = 0; //zerowanie sygnalu; x_k = 0 <=> syg[n] = 0
 }
 
//*****DEFINIOWANIE SYGNALU******************************************************************
 for(n = 0; n < len; n++)
 {
  if((n > 136) && (n < 148)) {syg[n]=250;} // dla danych przedzialow w zadaniu ustawiam sygnal 250
  if((n > 175) && (n < 211)) {syg[n]=250;} // f(x_k)=250 <=> syg[n]=250; k <=> n
 }
   
//******TRANSFORMATA FOURIERA****************************************************************
 for(n = 0; n < len; n++)
 {
  p->DFT[n].Re = 0.0;
  p->DFT[n].Im = 0.0;
  for(k = 0; k < len; k++)
  {
      p->DFT[n].Re +=  syg[k]*cos(2*M_PI*k*n/len);
      p->DFT[n].Im += -syg[k]*sin(2*M_PI*k*n/len);
  }
 } 
 
//******FILTR DOLNOPRZEPUSTOWY***************************************************************
/*
Amplituda jest rowna 250. 4% amplitudy jest rowne 10.

Po kilku probach dla roznych wartosci indeksu 
filtru dolnoprzepustowego (200,50,25,21,20) okazuje sie, 
ze oscylacje nie przekraczaja wartosci 10 
przy indeksie filtru 21 (i mniejszych).

Najwieksze oscylacje wystepuja przy pierwszym, krotszym sygnale, 
gdzie wahaja sie mniej wiecej w przedziale <240 data-blogger-escaped-259="">

(Arkusz FILTR 21)
*/
 for(n = 0; n < 22; n++) 
    {
     p->DFT[512-n].Re = 0.0;
     p->DFT[512+n].Re = 0.0;
     p->DFT[512-n].Im = 0.0;
     p->DFT[512+n].Im = 0.0;
    }
 
//******ODWROTNA TRNSFORAMTA FOURIERA********************************************************
 for(n = 0; n < len; n++)
    {
     p->iDFT[n].Re = 0.0;
     p->iDFT[n].Im = 0.0;
     
     for(k = 0; k < len; k++)
     {
         p->iDFT[n].Re += (p->DFT[k].Re*cos(2*M_PI*k*n/len) - p->DFT[k].Im*sin(2*M_PI*k*n/len));
         p->iDFT[n].Im += (p->DFT[k].Re*sin(2*M_PI*k*n/len) + p->DFT[k].Im*cos(2*M_PI*k*n/len));
        }
    }
    
//***********************************************************************************************
 
 pFile = fopen ("sygnal.txt", "wt");
    fseek(pFile, 0, SEEK_SET);
    
 for(i = 0; i < len; i++)
    {
     printf (    "%5d %4d %12.7f %12.7f %12.7f %12.7f\n", i, syg[i],  p->DFT[i].Re,  p->DFT[i].Im,  p->iDFT[i].Re,      p->iDFT[i].Im);  
     fprintf(pFile, "%5d %4d %12.7f %12.7f %12.7f %12.7f\n", i, syg[i],  p->DFT[i].Re,  p->DFT[i].Im,  p->iDFT[i].Re/(double)len, p->iDFT[i].Im/(double)len);
    }
    
    free(p);
    fclose(pFile);
 system("PAUSE"); 
 return 0;
}