Matlab Uygulamaları

Çok Noktanın Çok Noktaya Mesafe Hesabının Hızlı Yapılması, 15 Temmuz 2010
xy düzleminde iki grup noktamız olsun, ilk grupta N tane ikinci grupta M tane nokta olsun, bize lazım olan bu iki grup arasındaki noktaların herhangi ikisi arasındaki mesafe.  Örn: Birinci grubun t. noktasıyla ikinci grubun k. noktası arasındaki mesafe. Bunu içiçe iki for döngüsüyle yapmak mümkün ancak Matlab in vektörel işlem kabiliyeti sayesinde çok daha hızlı yapabiliyoruz. for döngüleri Matlab de yazılan kodların en zaman kaybettirici kısımları oluyor.

Birinci grubun x ve y lerine:
x1ler, y1ler diyelim.
İkinci grubun x ve y lerine:
x2ler, y2ler diyelim.

Bu dört diziyi kullanarak “mesafe” adında bir matris elde edeceğiz, bu matriste mesafe(t,k) birinci gruptan t. eleman ile ikinci gruptan k. eleman arasındaki mesafeyi döndürecek.
Kodumuz şöyle:
x1ler=[1 4 5 8 9 2 4];
y1ler=[2 3 7 2 6 2 1];
x2ler=[2 4 5 5 9 2 4 9 11 4 5];
y2ler=[3 2 5 6 4 1 8 0 1 3 1];
değerlerini veriyoruz, gerisi hesaplanıyor:
z1=x1ler+i*y1ler;
z2=x2ler+i*y2ler;

m1=repmat((1:length(z1)).’,1,length(z2));
m2=repmat(1:length(z2),length(z1),1);
mesafe=abs(z1(m1)-z2(m2));

Mesela yukarıdaki örnekte birinci gruptaki noktalar:
(1,2), (4,3), (5,7), (8,2), (9,6), (2,2), (4,1)
ikinci gruptaki noktalar:
(2,3), (4,2), (5,5), (5,6), (9,4), (2,1), (4,8), (9,0), (11,1), (4,3), (5,1)

mesafe(2,3) birinci grubun 2. elemanı ile ikinci grubun 3. elemanı arasındaki mesafeyi, yani: kok(1^2+2^2)=2.2361 verir, bizim bulduğumuz:
mesafe(2,3)
ans =
2.2361

Bu kodu tek grup noktalar için de kullanabiliyoruz, x1ler ve y1ler; x2ler ve y2ler ile aynı yapacağız.
Örn noktalarımız:
(2,3), (4,7), (6,2), (1,1) olsun.
x1ler=[2 4 6 1];
y1ler=[3 7 2 1];
x2ler=x1ler;
y2ler=y1ler;

>> z1=x1ler+i*y1ler;
z2=x2ler+i*y2ler;

m1=repmat((1:length(z1)).’,1,length(z2));
m2=repmat(1:length(z2),length(z1),1);
mesafe=abs(z1(m1)-z2(m2));
İlk noktanın kendisine uzaklığı:
>> mesafe(1,1)
ans =
0
1. notkanın 2.ye uzaklığı:
>> mesafe(1,2)
ans =
4.4721

N hücreli bir sistem düşünelim, bu N hücre birbirini aralarındaki uzaklığın bir fonksiyonu olarak etkiliyorsa yukarıdaki mesafe matrisi sayesinde NxN lik matris denklem sisteminin kurulması çok hızlanıyor. Benim kullandığım yer:
ros=abs(h_k(xindis)-h_k(yindis));
C=(1i*pi*k*a/2)*besselj(1,k*a)*besselh(0,2,k*ros).*repmat(contrast,length(h_k),1);
Hız >100/1 arttı.
İlgilisi olabilecek kesime sunarım.

Matlab FFT, 1 Haziran 2010
Matlab de FFT sonucunun işaretimizin gerçek frekans spektrumu gibi çıkması için birkaç düzenleme yapmamız gerekiyor. Matlab de FFT den düzgün sonuç alabilmemiz için Nyquist kriterinin en az 4-5 katı örnek almamız gerekiyor,  örnek sayısını bilgisayarınızı zorlamayacak kadar yüksek değerler seçebilirsiniz. Sinyalin gözlem süresini de uzun tutunca(en uzun periyodun 4-5 katı) FFT sonuçları daha sağlıklı oluyor.
Frkenas eksenimizi -1/(2*Ts) ile 1/(2*Ts) arasına ölçeklemeliyiz. Ts imiz örnekleme periyodu.  İnternetteki örneklerden de faydalanarak Matlab de FFT için küçük bir fonksiyon yazdım, adı pfft(proper fft):

function [] = pfft(t,x)
%proper fft, fatiherdem.net
%t, zaman vektörü
%x, örnek vektörü

T = t(2)-t(1);
fs = 1/T;
number_of_samples=length(x);
f = fs/2*linspace(-1,1,number_of_samples); 

fftSignal = fft(x);
fftSignal = fftshift(fftSignal);

subplot(1,2,1);
plot(t,x);
ylim([-2 2]);
xlabel('Time (t)');
ylabel('Magnitude');
subplot(1,2,2);
plot(f,abs(fftSignal));
%xlim([-1e6 1e6]); %xlim, x ekseni deger araligini beklediginiz frekans bandina gore
%secmelisiniz
%xlim([-1e4 1e4]);
xlim([0 1e4]);
xlabel('Frequency (Hz)');
ylabel('Magnitude');

end

İyileştirmeler(18.03.2014):
Genlik bilgisinin doğru çıması için:
plot(f,(2/length(t))*abs(fftSignal));
y eksenini ölçeklemek için:
var=max((2/length(t))*abs(fftSignal))*2;
ylim([0 var]);
yazılabilir.
Bazı örnekler:
1khz lik bir sinus düşünelim:

Periyodu 1ms oldu. 5 periyot boyunca gözleyelim ve örnekleme sıklığımız işaret frekansının 10 katı yani 10khz olsun, yani 0.0001 saniyede bir örnek alalım:
t=0:0.0001:0.005;
f=sin(2*pi*1000*t);
pfft(t,f);
fft1

10khz lik bir sinus ile 1MHz lik başka bir sinus işaretin çarpımına bakalım(DSB-SC AM modülasyonu):
Bu aslında DSB-SC AM(double side b. suppressed carrier AM) işareti oluyor, yani çift yanbantlı taşıyıcısı bastırılmış genlik modülasyonlu işaret, görüleceği üzere işaretimizde çift yanbant var ve taşıyıcısı görünmüyor.
En yüksek frekans 1MHz+10khz, 0.1us ile örnek alalım. En uzun periyotlu(0.1 ms) 10khz lik bileşen olduğundan 0.5ms boyunca sinyali gözleyelim, çünkü 10khz lik değişimi de yeterince gözlememiz gerekiyor:
t=0:1e-7:5e-4;
f=sin(2*pi*1e4*t).*sin(2*pi*1e6*t);
pfft(t,f);

pfft deki xlim sınırlarını:
xlim([0.9e6 1.1e6]);
olarak seçersek:
fft2

Tek yan bantlı(SSB) genlik modülasyonu:
1 Hz lik bir mesaj işaretini-tek ton-, 8 Hz lik taşıyıcı ile taşıyalım:
>> t=linspace(0,10,1000);
>> ton1=cos(2*pi*1*t);
>> ton2=sin(2*pi*1*t);
>> c1=cos(2*pi*8*t);
>> c2=sin(2*pi*8*t);
>> tyb_ust=ton1.*c1-ton2.*c2;
>> pfft(t,tyb_ust);
ssb_fft

Alt yan bant için tyb_alt=ton1.*c1+ton2.*c2; yazılması gerekir.

Testere işareti (sawtooth) ni fourier serisi toplamı olarak yazalım ve FFT sine bakalım:

f=1;
t=linspace(0,10,20000);
saw=0;
for n=1:999
saw=saw+(-1)^(n-1)*(2./(n*pi))*sin((2*pi*f)*n*t);
end
pfft(t,saw);

saw

 

Periyodik olmayan bir örnek, Darbe işareti:

t=-100:0.1:100;
y=zeros(1,length(t));
orta=round(length(y)/2);
y(orta-50:orta+50)=1;
pfft(t,y);
kare2

daha dar bir darbe:
t=-100:0.1:100;
y=zeros(1,length(t));
orta=round(length(y)/2);
y(orta-10:orta+10)=1;
pfft(t,y);
kare4
son iki örnekten zaman bölgesindeki sıkışmanın frekanstaki yayılmaya karşı düştüğünü görebiliyoruz.


Bu pfft fonksiyonunu ben emin olduğum spektrumları bir de FFT ile görmek için kullanıyorum, FFT sonucunun doğruluğu; sinyalin gözlem süresi, örnekleme sıklığı gibi parametrelere göre değişiyor. Sinyalimizi; sinyal içinde en yüksek frekansın 5-10 katı ile örnekleyip en uzun periyotlu sinyalin periyodunun 5-10 katı kadar gözlediğimizde sonucun doğruluğunu kuvvetlendirmiş oluyoruz. Henüz karşılaşmamış olsam da hatalı sonuçlar bulma riski bulunduğunu hatırlatayım.  İyi çalışmalar.

Tags :

2 thoughts on “Matlab Uygulamaları”

Bir yanıt yazın

E-posta adresiniz yayınlanmayacak. Gerekli alanlar * ile işaretlenmişlerdir