Cordic Algoritması Hakkında

[alicavuslu]

Açıkçası Hocam sayı formatı analizlerinin belirlenmesi işlemini VHDL yerine Matlab'da yapıyorum. Ama VHDL kodlamalarını yaparken generic olarak yazıyoruzki değişikliklerimiz kolay olsun.

FFT, Filtre, matris işlemleri vb uygulamalarda sayıların duyarlılıklarını Matlab ile ayarlıyorum. Donanıma yakın olması için fix komutu ile sayı yuvarlama işlemlerini yapıyorum ki bunu özellikle belirttim. round komutu ile tam olarak donanımsal işlem karşılığı olmuyor. Daha sonra floatin point formatında sonuçlarla karşılaştırıp ortalama karesel hata değerine bakıyoruz. Burada elde edilen hata oranı doğrultusunda sayı formatına karar verip generip parametresinde ilgili yerine yazıyoruz.

Cordic için böyle bir Matlab kodu olmadığından direk  simülasyon sonuçlarını verdim

[z]

Bende bu amaçla visual dillerle test programı yazıyorum.

real yada extended vs tanımlı değişkenlere atadığım sayıları 2^n ile çarpıp ardından round ile yuvarlıyorum.

Algoritmayı bu int değerlerle integer aritmetik yaparak işletiyorum. En son aşamada 2^-n ile sonucu elde edip bir de real değerlerle hesapladığım değer arasındaki farka bakıyorum.

N için 8, 16, 24, 32, 40, 48, 56, 64 ..... gibi değişik değerler için aradaki fark seçmiş olduğum bit sayısının yeterli olup olmadığını deneysel olarak elde etmemi sağlıyor.

Fakat bu çözümün matematiksel hiç bir özelliği yok.

Zaten sayısal sinyal işlemenin en can alıcı noktası burası. Burası yanlış kurulursa  inşa edilen yapı kesin çökecektir.

[alicavuslu]

Bu iş için açıkçası matematiksel bir yöntem var mı bilmiyorum. Ama daha önce hiç karşılaşmadım. Çünkü her sistem için farklı hata oranları sonuca olumsuz etkileyebilmekte. Bu yüzden sistem modelinin de haricen çıkarılıp analizler yapılması gerekebilir.

Mesela filtre çıkışına fft uygulayacaksanız mormalde binde birlik hata oranı filtre çıkışını çok etkilemez iken fft sonuçlarında hatalara neden olabiliyor. Dediğinzi gibi bunlarda yanlış yapıldığında yaptığınız işin bir önemi kalmıyor.

[z]

Bu bit sayısına deneysel yollardan ulaşıldığını düşünmüyorum kesinlikle bunun bir matematiği vardır.

Matematikçilerle anlaşamıyorum çünkü derdimi anlatamıyorum. Kendilerini vermiyorlar can kulağıyla sorumu dinlemiyorlar.



------------

1/(karekök(1+tan^2(x)))=cos(x) oluyor.

[alicavuslu]

Hocam analiziniz için teşekkür ederim.

Daha önce matematiksel yaklaşım olacağını hiç düşünmemiştim. Bit sayısının hesabı ile alakallı literatür taraması yapıp edindiğim bilgileri sizinle paylaşırım.

[z]

FIR sistemlerde bu belki çok fazla sorun değil ama IIR sistemlerde olmazsa olmaz.

Cevap bulursan sevinirim.

Aslında konuya bayağı kafa yormuşluğum var ama sonuca ulaşamadım. Kağıt kalemi kolay fırlatan tiplerdenim.


mesaj birleştirme:: 10 Eylül 2015, 17:27:09
Cordic algoritmasında eksik bir kısım kalmasın



Bu işlemlerde

Xn+1= Xn - Yn * Tan(DeltaTeta)
Yn+1= Yn + Xn * Tan(DeltaTeta)

Burada

Xn+1=sin Teta
Yn+1=cos Teta

Diyemiyoruz.

Diyebilmemiz için Bu değerleri 1/Xn ile çarpmamız gerekiyor.

Neden çarpmamız gerektiğini yukarıdaki işlemlerden görebilirsiniz.

Yn+1/Xn+1 oranında 1/Xn leri sadeleştirmiştik. Bu sadeleştirme tanjant değerini değiştirmez. Ancak bu pay ve payda daki değerlere sin ve cos dememize engel olur.

Her iterasyonda bir önceki iterasyonda sadeleştirilen 1/xn leri Xn  ile çarpmayı gerektirdiğinden

Cos(45) * cos(26.6) * cos(14) * cos(7.1) .... = 0.607252935 ye gidiyor.

yani cos(Atan(1))*cos(Atan(1/2))*cos(Atan(1/4))*cos(Atan(1/8))*........

[z]

Aliçavuşlu bir gün buluşalım sana kebap ısmarlıyayım.

Cordic algoritması

[JOKERAS]

z ve aliçavuşlu ustalar,

Hazır lafı geçmişken Şu filitreleri incelesek Fır,Iır,Ma vs.
Kaç çeşit filitre var,ne işe yararlar,nerelerde kullanlırlar vs.

Mesela bir ara CVD(Capasitive Voltage Divide) okumasıyla Capasitif matrix buton yapmak istemiştim.
Yaptım aslında ama tam istediğim gibi olmadı.
Bir kaç çeşit filitre algoritması yaptım yemedi,haddi leen dedi:)
Bu Adc okumalarında veya frekans değişim algılamalarında gürültü en büyük problem.

Şu filitre konusunu matematiksel olarak değilde C dili üzerinden tane tane mantığını anlatan biri olsa
Kebap ne kelime,Bursa İskender,buz gibi kola üstünede kaymaklı tel kadayıfı ısmarlarım:)
Varmısınız?Vaaarım diyooor:)

[z]

Digital filitreler

y(n)=K0x(n)+K1x(n-1)+K2x(n-2)+K3x(n-3)+.... şeklinde yazılırsa FIR filitre diye anılr.

Çıkış sadece girişin o anki ve geçmiş değerlerine bağlıdır.

Yok eğer

y(n)=K0x(n)+K1x(n-1)+K2x(n-2)+K3x(n-3)+.... + A0y(n)+A1y(n-1)+A2y(n-2)+A3y(n-3)+....şeklinde yazılırsa IIR filitre diye anılr.

Burada x(n) giriş sinyalinden alınan n. sample değeri anlamına gelirken y(n) n. sampledaki çıkış değeridir.

Digital filitre denen olay bundan ibaret.

Peki An ve Kn değerleri nasıl hesaplanacak.

Ben bu kısmı beceremiyorum. Bu yüzden Analog filtre tasarlayıp ardından bunu sayısallaştırıyorum.

n. dereceden bir filitre tasarımında, Matlab ya da özel filitre tasarım programları kullanmadan katsayıların hesaplanmasını anlatacak olan varsa ben de dinlemek isterim.

[JOKERAS]

Eyvallah z usta sağolasın.

Şimdi elimde sinüs 50Hz veya 100Hz daha yüksek frekanslardada olabilir böyle  bir sinyal var,bu sinyal temiz değil.
Sinyalin üzerinde bir sürü gürültü var,bu gürültüleri yok edip temiz bir sinüs elde etmek istiyorum.

Veya elimde 100KHz kare dalga bir sinyal var,bu sinyalin pulse'lerinde gürültü var ayrıca bu sinyal çevresel etkilerle
genel olarakta rastgele bir şekilde dalgalanıyor bunu nasıl filitre ederiz?

Yada Analog Audio sinyalide olabilir bunu dijitale dönüştürmek vs.


C dilinde örnek olursa süper olur.

[fatih6761]

Digital filitreler

y(n)=K0x(n)+K1x(n-1)+K2x(n-2)+K3x(n-3)+.... şeklinde yazılırsa FIR filitre diye anılr.

Çıkış sadece girişin o anki ve geçmiş değerlerine bağlıdır.

Yok eğer

y(n)=K0x(n)+K1x(n-1)+K2x(n-2)+K3x(n-3)+.... + A0y(n)+A1y(n-1)+A2y(n-2)+A3y(n-3)+....şeklinde yazılırsa IIR filitre diye anılr.

Burada x(n) giriş sinyalinden alınan n. sample değeri anlamına gelirken y(n) n. sampledaki çıkış değeridir.

Digital filitre denen olay bundan ibaret.

Peki An ve Kn değerleri nasıl hesaplanacak.

Ben bu kısmı beceremiyorum. Bu yüzden Analog filtre tasarlayıp ardından bunu sayısallaştırıyorum.

n. dereceden bir filitre tasarımında, Matlab ya da özel filitre tasarım programları kullanmadan katsayıların hesaplanmasını anlatacak olan varsa ben de dinlemek isterim.

Hocam FIR tasarımında özel bir şey yok diye biliyorum ben. Hatırladığım kadarıyla şöyle:

Bir filtre tasarlarken aslında freq. domain'de belirli bir alanı korurken diğer kısımları atmak istiyoruz. Bunu istediğimiz kısımları 1 ile diğerlerini 0 ile çarparak elde edebiliriz.
Freq. domainde bu çarpım fonksiyonunu plot ettiğimizde çift taraflı bir unit-step fonksiyonu oluyor. Yani parçalı fonksiyon olarak frekans;
-w0 dan küçükse 0
-w0 ile w0 arasında ise 1
w0 dan büyükse 0 dersek 0->w0 aralığını geçiren bir low pass filtre elde etmiş oluruz.
Bu fonksiyona ters fourier operatörünü uygularsak:
h(t) = integral(1*e^(-jwt)dw, from -inf, to +inf) olur. 0 olan kısımları integrale hiç dahil etmeye gerek yok.
İntegrali çözdüğümüzde;
h(t) = 2*sin(w0 * t) / t elde ederiz.

Freq. domainde (s-domainde olduğu gibi) çarpma işlemi time-domainde konvolüsyona karşılık gelir.
Time domain giriş sinyalimizi yani x(t) yi h(t) ile konvolüsyon işlemine tabii tutarsak sinyali teorik olarak filtrelemiş oluruz.
Peki discrete time da konvolüsyon nasıl olacak?
Bunun için filtre derecesini seçmemiz gerek. Mesela 7. dereceden olsun.
Bu 7 bize t için -3 ten 3 e kadar değerleri verir. Ayrık katsayılar
k[0] = h(-3)
k[1] = h(-2)
...
k[6] = h(3)
olacak şekilde 7 tane katsayı verir.
Giriş verimizin en güncel 7 değeri x[n] ... x[n-6] olsun.
Çıkış değerimiz konvolüsyona göre;
y[n] = x[n] * k[0] + x[n-1] * k[1] + x[n-2]*k[2] + ... + x[n-6]*k[6] olur.
Tabi burada continuous olan h fonksiyonunu discrete k'ya çevirirken oluşan veri kaybı filtrenin kazancında frekansa bağlı istenmeyen jump'lar oluşturur. Bunları azaltmak için windowing olayı uygulanır.
k değerleri hesaplanırken bir window functionın (triangle, hamming, kaiser vs..)  değerleriyle çarpılarak daha yumuşak geçişli bir katsayı listesi oluşturulur.

Bildiğim kadarıyla fir bu şekilde hocam. IIR yi hatırlamıyorum tam.

[z]

Şimdi anlayamadığım iki konu var.

1) Sen universiteye bu sene başlamış birisi olarak bunları nasıl bilebilirsin?

2)

Bunun için filtre derecesini seçmemiz gerek. Mesela 7. dereceden olsun.
Bu 7 bize t için -3 ten 3 e kadar değerleri verir. Ayrık katsayılar
k[0] = h(-3)
k[1] = h(-2)
...
k[6] = h(3)
olacak şekilde 7 tane katsayı verir.
Giriş verimizin en güncel 7 değeri x[n] ... x[n-6] olsun.
Çıkış değerimiz konvolüsyona göre;
y[n] = x[n] * k[0] + x[n-1] * k[1] + x[n-2]*k[2] + ... + x[n-6]*k[6] olur.

Bu kısmı biraz daha açabilirmisin?

[fatih6761]

Şimdi anlayamadığım iki konu var.

1) Sen universiteye bu sene başlamış birisi olarak bunları nasıl bilebilirsin?

2)
Bu kısmı biraz daha açabilirmisin?

1) özel ilgi diyelim hocam
2) Hocam asıl sorunuz filtre derecesiyle ilgili anladığım kadarıyla. Orada 0'ı ortaya almak için yapılıyor.
Aşağıda konvolüsyonu açıkladım, aslında sonsuz elemanla çarpma yapılması gerektiğinden ne kadar çok eleman kullanırsanız
o kadar iyi ama sonlu bir sayı olacak şekilde N. dereceden filtrede N (impulse response adım sayısına göre N) bir miktar hata olacak.
sinx/x fonksiyonu y eksenine göre simetrik olduğundan ve x=0 da en büyük değere sahip olduğundan göreli hatayı minimize etmek için 0'ı ortaya alacak şekilde soldan ve sağdan (N-1)/2 eleman alınıyor.
Yoksa 0..N veya -N..0 da alsak olur ama hata çok büyük çıkar. Burada hatadan kasıt freq. domeninde olduğundan filtre işlevini yitirecektir.

================================================================

Konvolüsyon teoremi: https://en.wikipedia.org/wiki/Convolution_theorem
Özellikle bizi ilgilendiren şu denklik:

Burada diyor ki;
"Frekans domeninde iki fonksiyonun çarpımı = Zaman domeninde o iki fonksiyonun konvolüsyonu"
Biz ne yapmak istiyoruz? Frekans domeninde bir fonksiyona belirli katsayılar uygulamak istiyoruz.
Freq. domeninde iki fonksiyonu çarpmak için normalde yöntem şu şekilde:
1) freq. domeninde katsayıları hesapla.
2) sinyali t. d. den freq. d. ye çevir.
3) fonksiyonları çarp
4) sonucu f.d. den t.d. ye geri dönüştür. Bu teoremi bulan da demişki bu kadar ara işleme gerek yok bu 4 basamak direkt time domain'de halledilebilir.
Teoremin ispatı konumuz kapsamında değil sanırım (wiki de verilmiş).

Peki konvolüsyon konvolüsyon diyoruz ne bu konvolüsyon? Convolution kelime anlamı olarak sarılma/sarmalanma demek yanılmıyorsam.
Genel tanım itibariyle cont. iki fonksiyonun konvolüsyonu (sarmalanması) şöyle hesaplanıyor:


Bizim durumumuzda sürekli olmayan ayrık (discrete) fonksiyonlar için bu integral Riemann toplamına dönüştürülebilir.

conv(f, g) =  sum ( f(n)*g(T-n) for n from 0 to N )

Olay bundan ibaret hocam.
Tabi kodlarken loop olarak yazıyoruz duruma göre. SSE AVX SIMD tarzı imkanlar varsa elle optimize ederiz ya da akıllı derleyiciler tree-loop-vectorizing denen bu işlemi bizim için yapar.

EDIT: Wiki'de de konvolüsyon güzel anlatılmış. Linki bulunsun:
https://en.wikipedia.org/wiki/Convolution

[z]

7 dereceden filitrenin katsayıları şunlarmı?

w0=1000 olsun.

h(t) = 2*sin(1000 * t) / t den

h(0)=2
h(1)=2sin(pi/3000)/3
h(2)=2sin(pi/3000)/2
h(3)=2sin(pi/3000)/1
diğer 3'ü de - işaretlileri.

[fatih6761]

7 dereceden filitrenin katsayıları şunlarmı?

w0=1000 olsun.

h(t) = 2*sin(1000 * t) / t den

h(0)=2
h(1)=2sin(pi/3000)/3
h(2)=2sin(pi/3000)/2
h(3)=2sin(pi/3000)/1
diğer 3'ü de - işaretlileri.

Direkt olarak değil hocam iki fark var:
Birincisi ifade tam olarak şu şekilde :
  h(t) = 2 * sinc(w0 * t) = 2 * ( sin (w0 * t) ) / t

Diğeri de frekans olarak normalized frequency kullanmamız gerekiyor. Yani w = 2*pi*f derken
f_norm = f_reel / f_sampling şeklinde. Yani sampling rate'e oranı şeklinde.
Sonra oradan w yı hesaplıyoruz.

Ekstra not olarak hocam hedeflenen gain'lere göre de katsayılar değişecektir logaritmik olarak.
Kağıt kalem olmadan kafamdakileri yazdığım için yanlış bi integrasyon bi katsayıyı atlama gibi hatalar olabilir.

He bir de bu arada şimdi farkettim filtrenin order'ı ile uzunluğu aynı şey değil. Yani 7 elemandan bahsediyorsak filter length= 7 oluyor.
Order = length - 1 yani 6.derece oluyor.

Evet kafadan integral almaya çalışınca hata oluyormuş gerçekten
w yerine fourier transorm için daha uygun olan 2*pi*f biçiminde ifade ederek integralini aldım: (wolframalpha)
http://www.wolframalpha.com/input/?i=integrate+e%5E%282*pi*i*f*t%29df+from+-f_0+to+f_0