Problem z FFT

0

Heja, mam problem z transforamta fouriera. Zaimpelementowalem algorytm radix 2 ale nie dziala on do konca dobrze. Tzn oglony charakter jest zachowany ale widac roznice w wartosciach. Ponizej przedstawiam program prawidlowy z matlaba i ten moj z C. Prosze o pomoc w znalezieniu bledu.

Matlab:

N=256;
x=0:N-1;
typBitReverse=2;
typFFT=2;
t=0:0.1:N/10-0.1;
x=sin(t*pi)+0.2*sin(t*pi*7);
xc=x;

a=1;

    for b=1:N-1

        if (b<a)
            T=x(a); x(a)=x(b); x(b)=T;
        end

        c=N/2;
        while (c<a)
            a=a-c; 
            c=c/2;
        end
        a=a+c;
    end
Moje progi w C:
// program realizuje test transforamaty fouriera i zapis do plikow wynikowych
#include <stdio.h>
#include <math.h>
#include "fourier.h"

#define N 256
#define LPI 3.14159265358979323846

int main(void)
{

    int i;
    double sinus[N],real[N],imag[N];

    FILE *plik_dane;
    FILE *plik_realis;
    FILE *plik_imag;    

//--------------------------------------------Utworzenie sygnalu sinusoidalnego    

    for(i=0;i<N;i++)
    {

           sinus[i]=sin(LPI*i*0.1)+0.2*(sin(LPI*i*7*0.1));

    }

   fourier(sinus,imag);

    //----------------------------------------zapis danych do pliku z tablic

    plik_realis=fopen("C:\\pliki/realis.txt", "w");

    for(i=0;i<N;i++)
       {
            fprintf(plik_realis,"%lf\n",sinus[i] );

       }

    fclose(plik_realis); 

    plik_imag=fopen("C:\\pliki/imag.txt", "w");

    for(i=0;i<N;i++)
       {
           fprintf (plik_imag,"%lf\n",imag[i] );

       }

    fclose(plik_imag); 

    return 0;
}

fourier.h

//realizacja fouriera dla 256 probek metoda radix-2

#define NN 256 //liczba probek
#include <math.h>

#define LPI 3.14159265358979323846

void fourier(double x[],double y[])
{
   int a=1,z,licznik;
   int b,c,d,e,g,L,M,mm,i;
   double T;
   double W_r,W_i,Wi_r,Wi_i,T_r,T_i;
   double licz[1024];

   licznik=0;

   FILE *wyn;

//-----------------------------------------przestawianie elementow

   for(b=1;b<NN;b++)
      {
        if (b<a)
            {
                T=x[a-1];
                x[a-1]=x[b-1];
                x[b-1]=T;
            }

        c=NN/2;

        while (c<a)
              { 
                 a-=c;
                 c/=2;
              }

        a+=c;
      } 

//--------------------------------------FFT

    for(e=0;e<NN;e++)
     { 
                    y[e]=0;
     }

    for(e=1;e<9;e++)
       {
           L = pow(2,e);
           M = L/2;
           Wi_r=1;
           Wi_i=0;

           W_r=cos(2*LPI/L);
           W_i=-(sin(2*LPI/L));

           for(mm=1;mm<M+1;mm++)           
              {
                  for(g=mm;g<NN+1;g+=L) 
                      {
                           d = g+M;

                           //T=x(d)*Wi

                           T_r = x[d-1]*Wi_r-y[d-1]*Wi_i;
                           T_i = x[d-1]*Wi_i+y[d-1]*Wi_r;

                            // x(d)=x(g)-T

                           x[d-1] = x[g-1]-T_r;
                           y[d-1] = y[g-1]-T_i;

                           //x(g)=x(g)+T

                           x[g-1] = x[g-1]+T_r;
                           y[g-1] = y[g-1]+T_i;

                      }   

              Wi_r=Wi_r*W_r-Wi_i*W_i;
              Wi_i=Wi_r*W_i+Wi_i*W_r;

              }
       }
}

Z gory dzieki za pomoc jest mi to naprawde potrzebne, pozdry

0

Z tego, co zauważyłam, przytoczony kod z Matlaba nie realizuje FFT, a jedynie odwrócenie bitowe próbek :-) W Matlabie FFT realizuje... funkcja fft(), a jeżeli chcesz zaimplementować ją ręcznie, to polecam książkę p. Tomasza Zielińskiego Cyfrowe przetwarzanie sygnałów. Od teorii do zastosowań.
Korzystając z tego kodu (w Matlabie) można łatwo napisać kod w C++, korzystając z klasy complex.

Pozdrawiam :-)

0

Fakt sory nie wkleilem reszty, oczywiscie mam ta ksiazke i lakgorytm wlasnie z tamtad pochodzi. Problem w tym ze ja to musze zrobic w C bo to do procka idzie a w Matlabie to kazdy potrafi fft wpisac, pozdrawiam:) A tu zalegly kod do matlaba:

for e =1 :log2(N)

    L = 2^e;
    M = 2^(e-1);
    Wi = 1;

    W = cos(2*pi/L)-j*sin(2*pi/L);

    for m = 1 : M
        for g = m : L : N

            f=f+1;
            d = g+M;
            T=x(d)*Wi;

            x(d) = x(g) - T;
            x(g) = x(g) + T;

           licznik(f)=d;

        end

        Wi=Wi*W;
    end

end
0

Pętla for(e=1;e<9;e++)
powinno być for(e=1;e<=9;e++)

0

Kilka uwagzamiast
L = pow(2,e);
M = L/2;
W_r=cos(2LPI/L);
W_i=-(sin(2
LPI/L));
napisałbym
L= 1<<e;
M= L/2;
W_r= cos(LPI/M);
W_i=-sin(LPI/M);

Potrzebne jest tylko log2(NN) wartości sinusa i cosinusa, możnaby pozbyć się obliczania tych funkcji.
Tablice z NN/2 wartości sin i cos pozwolą pozbyć się  zespolonego mnożenia wykonywanego w pętli po MM.

Zespolone mnożenie
c.r= a.rb.r - a.ib.i;
c.i= a.rb.i + a.ib.r
4-mnożenia, 2-dodawania,
można napisać inaczej,
t1 = a.r b.r;
t2 = a.i
b.i;
c.r = t1 - t2;
c.i = (a.r+a.i) * (b.r+b.i) - (t1+t2);
3-mnożenia, 5-dodawań, może być ważne gdy procesor to nie zmiennoprzecinkowy DSP


Współczesny kompilator powinien zamienić 
 for(mm=1; mm <M+1; mm++)    na
 for(mm=1; mm <= M; mm++) ale czy napewno?

Wszędzie przy indeksach do tablic masz minus jeden, pozbądź się tego.
0

Ok dzieki za rady. Juz dziala. Problem byl ze zwalilem zespolone mnozenia. A procesor to zmienno przecinkowy ad21061, w kazdym razie teraz jest idealnie jak w matlabie:)

0

Witam!! Wie ktoś może jak ten kod powyżej z C++ zapisać w Delphi? Pozdrawiam!!

0

wklej poprawny/uruchomiony kod to się zobaczy

0
mgr.Dobrowolski napisał(a)

Współczesny kompilator powinien zamienić
for(mm=1; mm <M+1; mm++) na
for(mm=1; mm <= M; mm++) ale czy na pewno?

moze, i bardzo czesto to robi, ale nie musi, wiec jesli mozna to warto samemu sie kolo tego zakrecic :) podobnie jest np. przy:

for(int i=0; i<kontener.size(); ++i)
//rob cos

pewnie jakies ponad 80% przypadkow nie zostanie to automatycznie zoptymalizowane do postaci

size_t size = kontener.size();
for(int i=0; i<size; ++i)
//rob cos

poneiwaz kompilator nie bedzie mial 100% pewnosci czy aby w //robcos sie wielkosc kontenera nie zmienia.. taka droba rzecz a duzo potrafi przyspieszyc. za 'starch czasow' za nie wywalenie takiego czegos poza petle sie oblewalo egzamin (->zasada wykluczania stalych poza iteracje), ale teraz na to juz zwracaja uwagi..

0

Hym , dlaczego pod Windows (XP) można poprawnie otworzyć plik
podając scieżki w formacie np :

 plik_realis=fopen("C:\\pliki/realis.txt", "w");

lub :

 plik_realis=fopen("C:/pliki/realis.txt", "w");

lub:

 plik_realis=fopen("C:\\pliki\\realis.txt", "w");

lub:

plik_realis=fopen("C://///pliki///realis.txt", "w");

itd:

plik_realis=fopen("C://////////////pliki/////////////realis.txt", "w");
  • jeszcze parę kombinacji .
    ?
    we wszystkich przypadkach plik zostaje poprawnie utworzony , wydawało mi się
    że jedynie przypadek 3 jest pod Win dopuszczalny ..
    Przepraszam że odbiegam od tematu .

1 użytkowników online, w tym zalogowanych: 0, gości: 1, botów: 0