FFT

[picusta]

Elime geçirdigim FFT algoritmasi.
Buyrun çözelim.
Cooley-Tukey  radical 2.
N boyunda bir FFT yapmak için N/2 boyunda 2 FFT yapiyoruz (nBlockSize).
Tek ve çift endeksli örnekleri bastan ayiriyoruz ve örnekleri yeniden siraliyoruz.
örnegin 8 örnekli bir FFT'de,FT  formülün dogrudan uygularsak 64 çarpma ve toplama yapiyoruz, FFT 'de ise 48. örnek sayisi artikça oran daha da artiyor:
Ekleme:
N örnekli bir Fourier dönüsümü için islem sayisi N² ile orantilidir, halbuki FFTde islem sayisi Nlog2 N ile orantilidir.
N = 1 milyon için, FFT kullanarak islem süresi 30 saniyedir, bu süre FT ile 2 haftayi bulur.
8 örnekli FFT islemi için asagidaki semayi ele alalim:

Sag tarafta ADC'den aldigimiz giris örneklerimiz var, x degiskeni bir sayi tablosu.
Sol tarafta ise X tablosu, yani Fourier dönüsümüne ugramis hali.
Giris örneklerinin arasinda D süresi varsa, tablonun i konumu iD zamanina denk geliyor. çikis sinyallerinin arasinda ise  i/ND frekans araligi vardir (N tablonun toplam uzunlugu).
FT'nin özellikleri arasinda, giris sinyali saf Reel ise sinyalin FT'si simetriktir.
Burada örneklerimiz adc'den çikmis saf reel oldugu için bu indirgeme yapilabilir (islem sayisi yariya iniyor).
1942'de Danielson ve Lanczos FT'nin baska güzel bir özelligini kesfettiler:
N boyunda bir FT yapmak için N/2 boyunda 2 FT yapip sonuçlari birlestirmek. Bu 2 FT'den birisi  N uzunlugudaki sinyalin tek endekslerinden olusmakta digeri ise çift endekslerden. Bu mantigi tekar edersek N/2 boyunda bir FT için 2 tane N/4 yapmak gerekir ... taa ki N=2'ye ulasana kadar. Peki 1 örnekli FFT ne verir? Girisi oldugu gibi çikisa verir.
Resimde görüldügü üzere algoritmanin basinda örnekleri öyle bir siraya yerlestiriyoruz ki, Fourier dönüsümünü özelliklerinden yararlanalim (asama, asama tek ve  çift endeksler) bunun için örneklerin tablodaki endeksini bitsel sirasini terseleyerek  degistiriyoruz.
1. Asamada nBlockSize=2, yani her blokta iki örnekli FFT yapmis gibi oluyoruz. Her asamada bu sayiyi 2 ile katliyoruz.
Referans kaynagi olarak "Numerical Recipe for C" kitabini kullanabiliriz, kitabin eski versyonu (1992) internette bulunuyor, ilgilendigimiz bölüm Fourier dönüsümü:
http://www.fizyka.umk.pl/nrbook/c12-0.pdf
http://www.fizyka.umk.pl/nrbook/c12-1.pdf
http://www.fizyka.umk.pl/nrbook/c12-2.pdf
http://www.fizyka.umk.pl/nrbook/c12-3.pdf
http://www.fizyka.umk.pl/nrbook/c12-4.pdf
http://www.fizyka.umk.pl/nrbook/c12-5.pdf
http://www.fizyka.umk.pl/nrbook/c12-6.pdf

static float WReal[MX_SAMPLE];
static float WImag[MX_SAMPLE];
void SetFFT(int N)
{
int   ii;
// W(n,k) katsayilarinin hazirlanmasi
for (ii=0; ii<N/2; ii++)
{
	WReal[ii] = cos(2*pi*ii/N);
	WImag[ii] = -sin(2*pi*ii/N);
}
}

void DoFFT(float *InReal, float *InImag, int N)
int ii,jj,					// Döngü degiskeni
	nBlockSize,				// Asamadaki Blok boyu
   wIndex,              // Frekans endeksi
	nWStep;					// Açisal (frekans) adimlama

float tempReal,
       tempImag;

// Giris sinyalini düzenle
for (ii=0; ii<N; ii++)
	{
	// yeni pozisyon :bitsel siranin terslenmis hali
	jj = reverseBits(ii,N);

	// iki sayinin yerini degistir
	if (ii<jj)
		{
		tempReal = InReal[ii];
		tempImag = InImag[ii];
		
		InReal[ii] = InReal[jj];
		InImag[ii] = InImag[jj];
		
      InReal[jj] = tempReal;
		InImag[jj] = tempImag;
		}          
	}
	
// FFT'nin hesaplanmasi
nWStep = N;
// Asama 
for (nBlockSize = 2 ; nBlockSize<=N ; nBlockSize*=2)
	{
	nWStep /= 2;
	for (i=0; i<N ; i+=nBlockSize)
		{
		wIndex = 0;
		k = i + n;
		j = k+ nBlockSize/2;
		
		wIndex += nWStep; // N/nBlockSize
		
		tempReal = WReal[wIndex] * InReal[j] - WImaag[wIndex]*InImag[j];
		tempImag = WReal[wIndex] * InImag[j] - WImaag[wIndex]*InReal[j];
		
		InReal[j] = InReal[k] - tempReal;
		InImag[j] = InImag[k] - tempImag;
		
		InReal[k] += tempReal;
      InImag[k] += tempImag;
		}

	}

[bymrz]

hocam fazla bişey anlamadık desek kızarmısınız...

hocam mesela bu uygulamada fiziksel olarak girdiler ve çıktılar nedir...
daha basit şekilde nasıl ifade edebilirsiniz...




Not: Ünv deki hocama çok kızıyorum şimdi fourier meselesini 5 dakikada bi grafik çizmişti anlatmış geşmişti, sizin zaten işinize yaramaz dercesine   :x

      Ve şunu öğrendim ki birisine birşey anlatmak ve anlamasınıda istiyosanız, ilk önce bunun ne işine yarayacağını nerelerde niçin kullanılacağını anlatmak lazım, yok rüzgar gibi gelip geçiyo...

Neyse fazla daldım:)

[kadirbas]

Selam picusta,

çözelimden kastınız, bu kodun aşama aşama ne yaptığını ya da formülsel olarak neyi ifade ettiğini mi öğrenmek, yoksa konuyu zaten bildiğinizi ve bu konuda bilgilendirici bir açıklama getireceğinizi mi söylüyorsunuz ??

[picusta]

Sabirsizlanmayin, aksam algoritmanin semasini scan edip koyarim, asama asama anlatilacak tabii.

[kadirbas]

Benimkisi sabırsızlıktan değil, yardım isteğindendi Kafanıza takılan bi noktada yardımım dokunabilir diye düşünmüştüm. Eskiden yazmıştım FFT kodları.

Yine de algoritma anlatımınızı bekliyor olacağım...

İyi çalışmalar.

[picusta]

FFT uygulamada baska sorunlar çikabiliyor.
Teorik olarak FT yapabilmek için sonsuz sayida örnek almak gerekir (-oo 'dan +oo'ya kadar entegre ediyoruz çünkü ), bunun imkansiz oldugu için örneklerimiz sinirli sayida. Bunu matematiksel olarak söyle açiklayabiliriz (Kisitli örnekli sinyal) = (Sonsuz örnekli sinyal) * (Pencere fonksyonu)
(Bu arada iki fonksyonun çarpiminin  FT'si konvolüsyon oldugunu belirtelim)
Bizim örneklemede pencere fonksyonumuz örnekleme süresince 1 olan ve diger yerlerde 0 olan bir fonksyon. Yani bir diktörtgen.

Solda pencere fonksyonumuzun zamana göre çizimi sagda ise bu fonksyonun FT dönüsümünün modülü.
ideal olarak pencere fonksyonumuzun FT'si bir dirak fonksyonu olmali ki konvolüsyon yapinca Kisitli örnekli sinyalimizin FT'si Sonsuz örnekli sinyalin FT'sine esit olsun.
Hamming pencere fonksyonun dönüsümü dirak'a yaklasiyor:

Hamming fonksyonun formülü :

Bu somut olarak demek oluyor ki, FFT yapmadan önce bütün örneklerimizi bu pencere fonksyonu ile çarparsak daha dogru sonuç aliriz.
Pencere fonksyonu olarak kullanilan baska bir fonksyon ise Hann:

[bymrz]

Sn kadirbas;

Bu konuyu bilmiyorum açıkcası, fakat öğrenmek istiyorum. Yani bu yüzden geliştirmeye yönelik bilgile değilde öncelikle öğrenmeye yönelik bilgilere ihtiyacım var..

Siz de yardımcı olursanız çok sevinirim...

Elinizde uğraştığınız veya yapmış oldugunuz uygulamalarda mevcutsa bizimle paylaşırsanız memun olurum..

[raltin]

hocam ellerine sağlık, bir süre önce araştırıp birşey anlamadığım bir konuydu . sayende bişeyler anlayabilicez, ama henüz çözemediğim bi kaç konu var :?
-saf reel derken sadece pozitif sayı anlamındamı?
-programı incelediğimde 2 dizi parametre alıyor, InReal ve InImag. bunlardan birisi tek endeksli birisi çift endekslilermidir?
-birde 16 bitlik integer için 1'e reverseBits uyguladığımda 32768 gibi bir değer çıkıyor,  bu doğrumudur (Algoritmayı delphi'ye aktarmaya çalışıyorumda delphide böyle bir değer çıkıyor), doğru ise elimizde en az 32768 örnekmi olması gerekir?

teşekkürler.
sürçi lisan eylediysem affola.

[picusta]

Haklisin Raltin, o fonksyon nitsel tersleme yapmiyor, bitlerin sirasini tersliyor.
örnegin 3 = 011b bitlerin sirasini degistirince : 110b = 6.
Ayni sekilde
0 = 000b  ->  000b = 0
1 = 001b ->   100b = 4
2 = 010b ->   010b = 2
3 = 011b ->   110b = 6
4 = 100b ->   001b = 1
5 = 101b ->   101b = 5
6 = 110b ->   110b = 3
7 = 111b ->   111b = 7
O fonksyonun prototipini söyle degistirebiliriz demek:

int reverseBits(int i, int N);

içini yazmak size kalmis. (N ile i sayisini kaç bitlik oldugunu bulup, 2 parçaya bölüp (low,high) sonra yerlerini degistirmek (öncesinden bir çikarma islemi gerçeklestirerek).)

[picusta]

-saf reel derken sadece pozitif sayı anlamındamı?
-programı incelediğimde 2 dizi parametre alıyor, InReal ve InImag. bunlardan birisi tek endeksli birisi çift endekslilermidir?

Basta sunu belirteyim, bu fonksyonda (DoFFT) InReal ve InImag tablolari, giris (semadaki soldaki dizi x sinyali) fonksyonun sonunda ise çikis (sagdaki dizi x sinyalinin dönüsümü X) degerlerini barindiriyor.

Saf reel derken durum söyle aslinda : FT yaparken giris sinyali komplex bir fonksyon olabilir. Ama bizim yapacagimiz uygulamada ADC'den bir sinyal ölçüp bunun spektrumunu elde etmek istiyoruz. Yani giriste komplex sayilar yok. sadece reel sayi dizisi var, örnegin :  145 + i*0, 178 + i*0, ...
InReal tablosu komplex sayinin reel kismini tutuyor. (145, 178, ...) InImag ise hayali kismi (0,0,...). Gerçekte uygulamamizda InImag tablosu (giris hayali kismi) hep 0'ile dolu olacak!
Kisacasi, verdigim kod sadece reel olan sayilar için tam optimize değil.
Optimize etmek için 2 yöntem var :
- N uzunlugundaki 2 reel fonksyonun FFT'sini hesaplamak için : InReal tablsouna 1. fonksyon, InImag tablosuna 2. fonksyon, daha sonra bu  2 fonksyonun FT'sini toplam FT'den ayirmak.
- 2N boyundaki reel fonksyonun  FFT'sini hesaplamak: InReel tablsouna çift endeksli örnekler, InImag tablosuna tek endeksli örnekler konarak, daha sonra çikis tablosunu ayirmak.

[einstain90]

Aslında çok basit bir örnek bulsak herkes daha rahat anlıcak.Mesela adc girişine atıyorum 100mhz sinyal verildiğinde bize bir değişkende  100 değerini veren bir örnek.Daha sonra bu frekansları ekrana yansıtmak çok kolay.

Mesela bir dizide  cikis[x] 

 cikis[0]=1  "100-299 mhz için " eğer cikis[0]=1(true) ise frekans aralıgı 100-299 dur
 cikis[1]=0  "300-399 mhz için "   .............
 cikis[2]=1   "400-499 mhz için " ..............
 cikis[3]=1   "500-599 mhz için " ..................
 cikis[4]=1   "600-399 mhz için "  .................

Şimdi sizin verdiğiniz programda çikiş

• hangisi? Giriş adc hangisi?Ama anladıgım kadarı ile biz şimdi girişe reel sayı veremiyoruz.

[picusta]

Girisler InReal ve InMag adreslerinden baslayan dizilerin içinde.
çikis ta ayni dizilerde yapiliyor.
girise saf reel sayi verebilirsin tabii, hiçbir engel yok.
Frekans ile tablonun endekslerini söyle bir formül bagliyor :
f = 2*pi*(örnekleme frekansi)  * Endeks / ( veri sayisi -1 )
Forumun programlama ve algoritma kismi oldugundan fazla kod vermek istemedim, sadece algoritma.
Projeler kisminda birisi spektrum cihazi yapsa söyle PIC+ grafik LCD üzerinde FFT hiç fena olmaz, karasimsek, kayan yazi projelerinden biktik

[einstain90]

peki ben adc olarak 125 okudum diyelim bunu inrealin hangi dizisine eşitlicem.
InReal[0]=125  InReal[1]=125........Bu böle gider.
Yazdıklarınız kafamı daha çok karıştırdı   Benim işlemci burda biraz zorlandı    

Yani işin mantıgını bir anlasam aninda lcdye dökücem.Ama giriş çıkışları bulamıyorum.Yani bana şu lazım adc neye eşitlicem sonucu hangi değişkende alıcam.Bunu anlasam olay bitti.Ama sanırım bunuda forumda kimse bulamıyo.
Bu konu daha öncede tartışılmıştı formuller verilmişti ama ondada bir proje yapamamıştık.

[picusta]

Basta N sayisini seçelim, spektrumdaki nokta sayisina esit olacak (256 diyelim).
örnekleme frekansi sabit. (5Khz farzedelim) demek oluyor ki program saniyede 5000 kere interrupta girecek.  
Her 256 örnekte bir, tablomuz dolmuş olacak. Iste bundan sonra fonksyonumuzu çagriyoruz ve 256 örnekli bir FFT hesaplamis oluyoruz.

[nthere]

Bu konuda burada kalmış.
Dökelim su FFTyi artık C kodlarına