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

Operacje na liczbach zespolonych – C/C++

// ©2015 Patryk Gradys
// OPERACJE NA LICZBACH ZESPOLONYCH

#include 
#include 
#include 

// STRUKTURA LICZBY ZESPOLONEJ ================================================== 

struct CMPL {
       double Re;
       double Im;
       };

// FUNKCJA DODAWANIA LICZB ZESPOLONYCH
CMPL C_SUM(CMPL z1, CMPL z2) {
         CMPL z3;
         z3.Re = z1.Re + z2.Re;
         z3.Im = z1.Im + z2.Im;
         return z3;
         }
         
// FUNKCJA ODEJMOWANIA LICZB ZESPOLONYCH
CMPL C_SUB(CMPL z1, CMPL z2) {
         CMPL z3;
         z3.Re = z1.Re - z2.Re;
         z3.Im = z1.Im - z2.Im;
         return z3;
         } 
           
// FUNKCJA MNOZENIA LICZB ZESPOLONYCH
CMPL C_MULT(CMPL z1, CMPL z2) {
         CMPL z3;
         z3.Re = z1.Re * z2.Re - z1.Im * z2.Im;
         z3.Im = z1.Re * z2.Im + z1.Im * z2.Re;
         return z3;
         }  
         
// FUNKCJA DZIELENIA LICZB ZESPOLONYCH
CMPL C_DIV(CMPL z1, CMPL z2) {
         CMPL z3;
         double w;
         w=z2.Re * z2.Re + z2.Im * z2.Im;
         
         if (w > 0) {       
             z3.Re = (z1.Re * z2.Re + z1.Im * z2.Im)/w;
             z3.Im = (z2.Re * z1.Im - z1.Re * z2.Im)/w;
             return z3;
         }
         else {
             printf("FUNCTION C_DIV ERROR: Wystapil blad podczas dzielenia \n");
             system("PAUSE");
             exit(10);
             }
         }                  
          
// =============================================================================
// FUNKCJA GLOWNA ==============================================================

int main () {
    
    CMPL z1, z2, z3, z4, z1_7, z2_9, z3_3, z4_5, licznik, mianownik, wynik;
    int i = 0;
    
    z1.Re = 1.0; // zmienna dla z1 (przypisywanie wartosci z zadania)
    z2.Re = 1.0;
    z3.Re = 2.0;
    z4.Re = 0.2;  
    
    z1_7.Re = 0.0; // zmienna dla z1^7 (wstepne zerowanie)
    z2_9.Re = 0.0;
    z3_3.Re = 0.0;
    z4_5.Re = 0.0;
    
    licznik.Re = 0.0; // zmienna dla licznika ulamka
    mianownik.Re = 0.0; // zmienna dla mianownika ulamka
    wynik.Re = 0.0; // zmienna dla wyniku
    
    z1.Im = 1.0; //przypisywanie wartosci z zadania
    z2.Im = -0.2;
    z3.Im = -1.0;
    z4.Im = -0.3;
    
    z1_7.Im = 0.0;  
    z2_9.Im = 0.0;
    z3_3.Im = 0.0;
    z4_5.Im = 0.0;
    
    licznik.Im = 0.0;
    mianownik.Im = 0.0;
    wynik.Im = 0.0;
    
    //obliczanie poteg liczb zespolonych
    
    // z1^7         
    z1_7 = z1;
    for(i = 1; i < 7; i++) {
          z1_7 = C_MULT(z1, z1_7);
          }   
    printf("z1^7:\t %8.4f \t %8.4f \n", z1_7.Re, z1_7.Im);           
    
    // z2^9    
    z2_9 = z2;
    for(i = 1; i < 9; i++) {
          z2_9 = C_MULT(z2, z2_9); 
          }   
    printf("z2^9:\t %8.4f \t %8.4f \n", z2_9.Re, z2_9.Im);     
          
    // z3^3          
    z3_3 = z3;
    for(i = 1; i < 3; i++) {
          z3_3 = C_MULT(z3, z3_3); 
          }  
    printf("z3^3:\t %8.4f \t %8.4f \n", z3_3.Re, z3_3.Im);     
       
    // z4^5      
    z4_5 = z4;
    for(i = 1; i < 5; i++) {
          z4_5 = C_MULT(z4, z4_5);     
          }  
    printf("z4^5:\t %8.4f \t %8.4f \n", z4_5.Re, z4_5.Im);  
          
    // mnozenie liczb z licznika
    licznik = C_MULT(z1_7, z2_9); 
    printf("licznik:\t %8.4f \t %8.4f \n", licznik.Re, licznik.Im);  
    
    // odejmowanie liczb z mianownika  
    mianownik = C_SUB(z3_3, z4_5);     
    printf("mianownik:\t %8.4f \t %8.4f \n", mianownik.Re, mianownik.Im); 
    
    //obliczanie ulamka (dzielenie licznika i mianownika)
    wynik = C_DIV(licznik, mianownik);
    printf("\nWYNIK:\t %8.4f \t %8.4f \n\n", wynik.Re, wynik.Im);   
    
    /*
    z3 = C_SUM(z1, z2); // wywolanie funkcji C_SUM
    printf("Dodawanie \n");
    printf("%7.2f %7.2f \n", z1.Re, z1.Im);
    printf("%7.2f %7.2f \n", z2.Re, z2.Im);
    printf("%7.2f %7.2f \n", z3.Re, z3.Im);
    
    z3 = C_SUB(z1, z2); // wywolanie funkcji C_SUB
    printf("Odejmowanie \n");
    printf("%7.2f %7.2f \n", z1.Re, z1.Im);
    printf("%7.2f %7.2f \n", z2.Re, z2.Im);
    printf("%7.2f %7.2f \n", z3.Re, z3.Im);
    
    z3 = C_MULT(z1, z2); // wywolanie funkcji C_MULT
    printf("Mnozenie \n");
    printf("%7.2f %7.2f \n", z1.Re, z1.Im);
    printf("%7.2f %7.2f \n", z2.Re, z2.Im);
    printf("%7.2f %7.2f \n", z3.Re, z3.Im);
    
    z3 = C_DIV(z1, z2); // wywolanie funkcji C_DIV
    printf("Dzielenie \n");
    printf("%7.2f %7.2f \n", z1.Re, z1.Im);
    printf("%7.2f %7.2f \n", z2.Re, z2.Im);
    printf("%7.2f %7.2f \n", z3.Re, z3.Im);
    
// POTEGOWANIE LICZB ZESPOLONYCH
    int i;
    CMPL z4 = z1;
    printf("Potegowanie \n"); 
    for(i = 0; i < 19; i++) {
          z4 = C_MULT(z1, z4);
          
    // Przy jakim indeksie modul potegi liczby zespolonej bedzie wiekszy niz 1024   
    printf("%3d %10.2f %10.2f %10.3f \n", i, z4.Re, z4.Im, sqrt(z4.Re*z4.Re + z4.Im*z4.Im));   
          }
    */  
    
    system("PAUSE");
    return 0;
}

Czytaj dalej...