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

Czytaj dalej...

Funkcje trygonometryczne na przedziale z zadanym rastrem – C/C++

// ©2015 Patryk Grądys
// Funkcje trygonometryczne na przedziale z zadanym rastrem

#include 
#include 
#include 

int main () {

    int i;
    double x, wynik;
    wynik = 0.0;
    double eps = 0.0000001;
    FILE *pFile;
    
    pFile = fopen ("prog2.txt", "wt");
    fseek(pFile, 0, SEEK_SET);

    printf("   i \t x \t wynik \n\n");         
        
    for(i=0; i<=100; i++) { //ustawiamy 100 wartosci
             x = 2.0*i/100.0-1.0; //przeskalowanie
             
             // bezpiecznik na wypadek dzielenia przez zero, 
             
             if(fabs(x) < eps) {                 
                  wynik = 999999.0; //wowczas wynik = 999999.0 (zeby byl latwo rozroznialny)
             }
             
             else {
             //sin^2(x) = (1-cos(2*x))/2 // z tozsamosci trygonometrycznej
             // Przeksztalcam wzor funkcji na jezyk C:
                              
                  wynik = ((1-cos(2.0*x))/2.0) / (sqrt(1.0-cos(2.0*x)));
             }
                 
             // wyswietlanie wynikow z przedzialu <-1 data-blogger-escaped-1="">
             if (i <= 100) { 
                 printf(          "%4d %8.4f %8.4f \n", i, x, wynik);  
                 fprintf(pFile,   "%4d %8.4f %8.4f \n", i, x, wynik);
             }  
             // co prawda w zadaniu jest mowa o odczycie tylko z dysku, 
             // ale wyswietlam pogladowo, zeby sprawdzic czy aby na pewno 
             // wlasciwie wartosci znajduja sie w pliku       
    }
             
    fclose( pFile);
    
    printf("\n Dane wyswietlane z pliku \n");
    
    // Wyswietlanie danych z pliku
    pFile = fopen ("prog2.txt", "rt");
    fseek(pFile, 0, SEEK_SET);
    
    int k;
    float y, fwynik;
    
    for(i = 0; i<=100; i++) { //wyswietlam 100 wartosci z pliku
          fscanf(pFile, "%d %f %f \n", &k, &y, &fwynik);
          printf("%4d %8.4f %8.4f \n", i, y, fwynik);      
    }
    
    fclose( pFile);
    
    system("PAUSE");
    return 0;
}

Czytaj dalej...