Matlab Sinyal Işleme

1 SİNYALLER VE SİSTEMLERİN MATLAB YARDIMIYLA BENZETİMİ Sinyal Üretimi Bu laboratuarda analog sinyaller ve sistemlerin sayısal bir ortamda benzetimini yapacağımız için örneklenmiş sinyaller üzerinde işlem yapacağız. Sürekli zamanlı bir sinyal olan X c ( t ) yi N boyutlu bir vektör ile örneklenmiş ve kuantalanmış (bölütlenmiş) olarak [X c(t 0) X c(t 1) X c(t N-1)] şeklinde ifade edilebilir. Burada sinyalin zaman aralığı t o ile tn 1 arasındadır. Örnekleme aralığı T s=ti+1-t i ifadesiyle belirtilmektedir. Örnekleme aralığı yeteri kadar büyük seçilmelidir ki MATLAB de sinyal sürekli zaman gibi görülebilsin. Sinyalin en büyük frekansının yaklaşık 10 katı bir değer (örnekleme frekansı) işlemler için yeterli olacaktır. Ancak sinyalin fazının çizdirilmesi işleminde alınacak değer katı olması yapılacak işlemin doğruluğunu artıracaktır. Örnekleme frekansı ile örnekleme zaman aralığı arasındaki bağıntı f = 1/ T dir. s s fs = 10 Hz lik bir sinüzoidal sinyal aşağıdaki gibi üretilir: %program ch2_1.m close all % Ekranda daha önce çizilmiş şekil varsa bu şekilleri kapatır. clear all % Daha önceden yapılmış bir işlem varsa hafızayı temizler. clc % Komut penceresi ekranını temizler. fm=10; % İşaretin frekansı 10 Hz fs=*fm; %Sinyalin örnekleme frekansı Hz; ts=1/fs; n=[0:(1/fs):1]; % Sinyal 0&#;dan 1 saniyeye kadar faz=0; %30 tsy=cos(2*pi*n*fm+faz); % İşaretimiz plot(n,tsy, &#;k&#; ); %işaretin zaman izgesinde çizimi title(&#;cosinus dalgasi&#;) xlabel(&#;saniye&#;); ylabel(&#;genlik&#;); 12

2 Şekil Kosinüs Dalgası Şekilden görüleceği üzere 1 saniyede 10 adet kosinüs dalgası vardır Fourier Dönüşümü Fourier dönüşüm yöntemi sinyalin içindeki bilgilerin elde edilebilmesi için, sinyallerin işlenmesinde kullanılan çok önemli bir yöntemdir. Bu bilgiler, Fourier dönüşümü ile MATLAB tarafından yeniden kullanılmaya uygun bir veri formatına çevrilir. Fourier dönüşümüyle bir sinyal, farklı genlik, frekans ve fazlarda kosinüs ve sinüs temel bileşenlerinin toplamı olarak ifade edilir. Her bileşenin frekans ve genliği ile birlikte tablolaşması, bilgisayarla verilerin işlenmesi sırasında kolaylık sağlar. + jw jwt X ( e ) = Xc ( t) e dt jw X( e ) = x( n) e jwn () () Denklem () Sürekli Zamanlı Fourier Dönüşümü, Denklem () ise Ayrıklı Zamanlı Fourier Dönüşümü nü göstermektedir. MATLAB de kullandığımız dönüşüm ise hem zamanda hem de frekansta ayrık olduğu için DFT ve IDFT kullanırız. N 1 jw 2 kn π X ( k) = x( n) e wk = k n= 0 N N 1 1 jwk n xn ( ) = X( ke ) N n= 0 () () Denklem () Ayrık Fourier Dönüşümü (DFT), Denklem () ise Ters Ayrık Fourier Dönüşümü (IDFT) dür. 13

3 Bu dönüşüm hesaplamaları maalesef çok masraflı hesaplamalardır. Hızlı Fourier dönüşümü tekniği, bir yandan hesaplamalar sürerken, bir yaklaşım olarak ilk elde edilen değerlerin kullanıma sunulduğu bir alternatif yazılım tekniğidir. %program ch2_2.m fm=10; % Isaretin frekansi fs=*fm; %Sinyalin ornekleme frekansi ts=1/fs; n=[0:(1/fs):1]; % Sinyal 0&#;dan 1 saniyeye kadar faz=0 %;%30 tsy=sin(2*pi*n*fm+faz); % isaretimiz % cos((pi/2)-a)=cos(a-(pi/2)) tsyf=fft(tsy)/length(tsy); % Sinyalin frekans izgesinde gösterilimi tsyfm=abs(tsyf); % Sinyalin fourier dönüşümü yapılınca karmaşık %faz temizle% for i=1:(fs+1); if abs(tsyf(i))<; tsyf(i)=0; tsyfa=angle(tsyf); SUBPLOT(2,1,1) plot(tsyfm); AXIS([ ]) title(&#;sekil 2a&#;); %text(,,&#;-pi +pi araligi için &#;) SUBPLOT(2,1,2) plot(tsyfa); AXIS([ ]) title(&#;sekil 2b&#;); xlabel(&#;hertz&#;); % Sinyalin frekansının bulunması [A,B]=max(tsyfm(1:(fs/2))); disp(&#;sinyalin frekansi&#;) disp(b-1) tsyfa(b) Şekil a. da sinyalin frekans cevabının mutlak değeri çizilmişken Şekil b. de faz cevabı çizilmiştir. Burada dikkat edilmesi gereken husus faz cevabı bulunurken örnekleme frekansının yeterince büyük seçilmesinin gerekliliğidir. 14

4 Şekil a) Sinyalin Frekans Cevabının Mutlak Değeri b) Sinyalin Faz Cevabı Sinyalin Fourier dönüşümünden sonra sıfıra yakın sayılar oluşmaktadır. Bu sayıların oluşumundan dolayı faz cevabı anlaşılır şekilde çıkmamıştır. Bu problemin çözülebilmesi için DFT işleminden sonra sıfıra yakın sayılar sıfırlanır. Bu işlem verilen örnekte faz temizleme ile kısmında gerçekleştirilmiştir Süzgeç Yapıları Süzgeçler yapılarına göre Sonlu Dürtü Yanıtlı (FIR, Finite Impluse Response) ve Sonsuz Dürtü Yanıtlı (IIR, Infinite Impluse Response) süzgeçler olmak üzere ikiye ayrılır. Ayrıca süzgeçler verdikleri frekans cevabına göre Alt Geçiren (LP, Low Pass), Üst Geçiren (HP, High Pass), Band Geçiren (BP, Band Pass), Band Bastıran Geçirmeyen (BR, Band eject), Tüm Geçiren (AP, All Pass) şeklinde ayrılmaktadır Sonsuz Dürtü Yanıtlı Süzgeç Yapıları Bu bölümde Butterworth süzgeç tasarımı verilecektir. Örnekte verilen tasarım alt geçiren süzgeç tasarımıdır ancak % li kısımlar kaldırılarak diğer tasarımların da nasıl yapılabileceği görülebilir. 15

5 %program ch2_5.m close all wg=[]; wd=[]; %wg=[] %wd=[] %wg=[ ]; %wd=[ ]; %wg=[ ]; %wd=[ ]; gddb=1; sddb=40; % Alt geçiren % Üst geçiren % Band geçiren % Band Durduran [N,Wn]=buttord(wg,wd,gddb,sddb); [B,A] = BUTTER(N,Wn); fs=; [H,W] = FREQZ(B,A,); Hg=20*log10(abs(H)); plot(w/pi,hg) grid on AXIS([ ]) xlabel(&#;w/pi&#;); ylabel(&#;kazanç,db&#;); title(&#;iir,buttordworth Alt Geçiren Süzgeç&#;) plot(abs(h)); grid on xlabel(&#;hz&#;); ylabel(&#;kazanç&#;); AXIS([ ]) for i=1:(length(h)); if abs(h(i))<; H(i)=0; Ha=angle(H); plot((ha/pi)*); xlabel(&#;hz&#;); ylabel(&#;faz&#;); grid on 16

6 Şekil Butterworth Alt Geçiren Süzgecin Kazanç Eğrisi (db) Şekil Butterworth Alt Geçiren Süzgecin Kazanç Eğrisi 17

7 Şekil Butterworth Alt Geçiren Süzgecin Faz Cevabı Ayrıca Chebyshev, Elliptic süzgeçleri de verilen MATLAB fonksiyonları ile gerçekleştirilebilir. buttord: Geçiş bandında R p (db) değerinden fazla olmayan ve durdurma bandında en az R s (db) değeri kadar güç yetirimini veren en düşük dereceli sayısal Butterworth süzgecin derecesini verir. W g ve W d geçiş ve durdurma bandının 0 ile 1 arasında normalize edilmiş köşe frekanslarını göstermektedir. Fonksiyonun çıkışı olan W n ise istenen özellikte süzgeç için gerekli olan doğal frekansı vermektedir. Alt Geçiren: W p = W s =; Üst Geçiren: W p = W s =; Band Geçiren: W p =[ ], W s =[,]; Band Durduran: W p =[ ], W s =[ ]; Butter: Butterworth sayısal ve analog süzgeç tasarımı N. dereceden alt geçiren süzgeç tasarlar ve N+1 uzunluğunda B (pay) ve A (payda) süzgeç katsayılarını verir. Katsayılar kaydırmalı yapı düşünülerek z in sıfırın kuvvetinden N+1. kuvvetine kadar gider. Ayrıca kesim frekansı W n < W n < arasındadır. Burada örnekleme hızının yarısını göstermektedir. Eğer W n iki bileşen oluşuyorsa W n = [W 1 W 2 ] 2N dereceli geçiş bandı W 1 < W < W 2 şeklinde olan süzgeç olur. Ayrıca üst geçiren süzgeç [B,A] = butter(n,wn,&#;high&#;) ile bant durduran süzgeç ise [B,A] = butter(n,wn,&#;stop&#;) ile tasarlanabilir. cheb1ord: Birinci çeşit Chebyshev süzgeç derecesi bulma cheby: Chebyshev birinci çeşit sayısal ve analog süzgeç tasarımı cheb2ord: İkinci çeşit Chebyshev süzgeç derecesi bulma cheby2: Chebyshev ikinci çeşit sayısal ve analog süzgeç tasarımı ellipord: Elliptic süzgeç derecesi bulma ellip: Elliptic veya Cauer sayısal ve analog süzgeç tasarımı 18

8 Sonlu Dürtü Yanıtlı Süzgeç Tasarımı FIR süzgeçlerin getirisi doğrusal faz cevabına sahip olmalarıdır. Ancak bu süzgeç yapılarında istenen frekans cevabını elde etmek için gerekli olan süzgeç uzunluğu oldukça fazladır. Dahası FIR tasarımında geçiş bandı ile durdurma bandı arası olan dönüşüm bandı IIR süzgeçler kadar keskin olmamasıdır. FIR süzgeçlerde çeşitli tasarım metotları vardır. Bunlar pencereleme, remez algoritması, en küçük kareler yöntemi gibi çeşitli yöntemlerdir. Bu bölümde bir pencereleme yöntemi ile yapılan süzgeç tasarımını vereceğiz. Bu tasarım varsayılan olarak hamming pencereleme yöntemini kullanmaktadır. %program ch2_6.m close all Wn=[]; % Alt geçiren N=; B = FIR1(N,Wn,&#;low&#;) fs=; [H,W] = FREQZ(B,1,); Hg=20*log10(abs(H)); plot(w/pi,hg) grid on AXIS([ ]) xlabel(&#;w/pi&#;); ylabel(&#;kazanç,db&#;); title(&#;sdc,hamming Pencereleme Yöntemiyle Alt Geçiren Süzgeç Tasarimi&#;) plot(abs(h)); grid on xlabel(&#;hz&#;); ylabel(&#;kazanç&#;); AXIS([ ]) for i=1:(length(h)); if abs(h(i))<; H(i)=0; Ha=angle(H); plot((ha/pi*)); xlabel(&#;hz&#;); ylabel(&#;faz&#;); grid on 19

9 Şekil FIR Alt Geçiren Süzgecin Kazanç eğrisi (db) Şekil FIR Alt Geçiren Süzgecin Kazanç Eğrisi 20

10 Şekil FIR Alt Geçiren Süzgecin faz cevabı Unwrap: Daha önce FIR süzgeç yapısı doğrusal bir faz cevabı verirken IIR yapılar bu cevabı veremez demiştik. Bu ifade çizdirilen faz grafiklerinde tam olarak görülememektedir. Bu yüzden unwrap denilen π den büyük atlama fazlarını 2π nin katlarına dolayan işlev kullanılır ve sürekli hali görüntülenebilir. Ha=angle(H); komutundan sonra Ha=unwrap(Ha); komutu kullanılarak yapılırsa sonlu ve sonsuz darbe cevaplı süzgeçler için aşağıdaki şekiller elde edilebilir. Buradan görülebileceği üzere sonsuz darbe cevaplı süzgecin faz eğrisi doğrusal iken bu eğri sonlu dürtü cevaplı süzgeç için doğrusala yakın ancak doğrusal değildir. Şekil Sonlu Dürtü Cevaplı Süzgeç İçin Faz Cevabı Düzenlenmiş Eğri 21

11 Şekil Sonsuz Dürtü cevaplı Süzgeç İçin Faz Cevabı Düzenlenmiş Eğri Sinyallerim Süzgeçlenmesi yt ( ) = funduszeue.info(2 π f t+ θ ) + funduszeue.info(2 π f t+ θ ) şeklinde verilen y() t 1 1c 1 2 2c sinyalini A = 3, A = 5, f = 10 Hz, f = Hz, θ = 30, θ = 0 değerleri için benzetimini yaparak daha sonradan yk() t = funduszeue.info(2 π f1 ct) sinyalini elde etmek için alt geçiren sonlu cevaplı süzgeçten geçirelim. Dikkat edilmesi gereken husus süzgeç tasarımı yaparken süzgeç tasarım kriterini 0 ile π arasında göz önüne alınmasıdır. Örneğin örnekleme frekansımız olsun ve kullanacağımız süzgeç bant geçiren olsun ve geçirme frekansları da ile Hz arasında olsun, bu durumda durdurma frekanslarını da ile seçelim. Süzgecimizin en büyük frekansı f s /2 olacak şekilde almamız gereken değerler basit bir oran orantı w = olarak bulunur. Burada önemli bir nokta bu ile w [ ] =, g d [ ] tasarımda bile seçtiğimiz değerlere karşılık gelen B değerinin çok küçük olmasıdır. Bu açıdan dar bantlı bir süzgeç tasarımının ve gerçekleştirilmesinin zor olduğu görülmektedir. Örneğin tasarımda B nin değeri direkt sıfır olarak alınırsa süzgeç çalışmaz. % Sinyal uretimi ve filtreleme islemi close all clear all hold on fs=; % pi= %Sinyalin ornekleme frekansi ts=1/fs; n=[0:(1/fs):1]; % Sinyal 0&#;dan 1 saniyeye kadar fm1=10; % fm2=; faz1=(pi/6); faz2=0; A1=4; A2=3; % Isaretin frekansi 22

12 tsy=a1*cos(2*pi*n*fm1+faz1)+a2*sin(2*pi*n*fm2+faz2); % isaretimiz %gy2=a2*sin(2*pi*n*fm2+faz2); süzgeçleme isleminden sonra istenen sinyal %Sekil 1: sinyalin zaman izgesinde çizimi plot(tsy); AXIS([ ]); xlabel(&#;örnek sayisi, Toplam saniye&#;); ylabel(&#;genlik Degeri,Volt&#;); % tsyf=fft(tsy)/length(tsy); tsyfm=abs(tsyf); %faz temizle% for i=1:(fs+1); if abs(tsyf(i))<; tsyf(i)=0; tsyfa=angle(tsyf); % % pi araliginda olan izgeyi -pi +pi araligina goturme islemi tsyfm=fftshift(tsyfm); tsyfa=fftshift(tsyfa); eks=[-fs/fs/2]; %Süzgeç Tasarimi(Burada süzgeçin frekans çiziminin %gösterimi için gerekli islemler yapilmaktadir.) wg=[ ]; % Band geçiren wd=[ ]; gddb=1; sddb=40; [N,Wn]=buttord(wg,wd,gddb,sddb); [B,A] = BUTTER(N,Wn); %B=0; [H,W] = FREQZ(B,A,fs/2+1); eh=flipud(h); H=[eH(1:fs/2);H]; % % subplot(2,1,1); plot(eks,abs(h));grid on AXIS([ ]) %faz temizle 23

13 for i=1:(fs+1); if abs(h(i))<; H(i)=0; subplot(2,1,2) plot(eks,angle(h)),grid on AXIS([ ]); xlabel(&#;hertz&#;); % %y(t) sinyali&#;nin cizimi SUBPLOT(2,1,1) plot(eks,tsyfm); hold on plot(eks,abs(h));grid on AXIS([ ]) hold on SUBPLOT(2,1,2) plot(eks,(tsyfa/pi)*); AXIS([ ]); hold on % plot(eks,angle(h));grid on %AXIS([ ]); xlabel(&#;hertz&#;); %Süzgeçleme islemi suz_tsy=filter(b,a,tsy); % suz_tsyf=fft(suz_tsy)/length(tsy); for i=1:(fs+1); if abs(suz_tsyf(i))<; suz_tsyf(i)=0; suz_tsyfa=angle(suz_tsyf); suz_tsyfm=abs(suz_tsyf); suz_tsyfm=fftshift(suz_tsyfm); suz_tsyfa=fftshift(suz_tsyfa); %Sekillerin çizdirimi SUBPLOT(2,1,1) plot(eks,suz_tsyfm);grid on AXIS([ ]); 24

14 SUBPLOT(2,1,2) plot(eks,(suz_tsyfa/pi)*); AXIS([ ]);grid on xlabel(&#;hertz&#;) plot(suz_tsy); gy2=a2*sin(2*pi*n*fm2+faz2); hold on plot(gy2,&#;r&#;); AXIS([ ]); xlabel(&#;örnek Sayisi, Toplam saniye&#;) Şekil y(t) işareti Şekil Kullanılan Bant Geçiren Süzgecin Kazanç ve Faz Cevabı 25

15 Şekil Y(F) İn Genlik Ve Faz Cevabı Üzerine Kullanılan Süzgecin Genlik Ve Faz Cevabının Gösterimi Şekil Süzgeçleme İşleminden Sonra İşaretin Frekans Cevabı 26

16 Şekil Süzgeçlenmiş Ve Gerçek İşaretin Gösterimi 27