Modul Pelatihan Matlab 2014

Modul

Pelatihan MATLAB Aplikasi Bidang Geofisika

Oleh : Dr. Andri Dian Nugraha Rexha Verdhora Ry, ST. Waskito Pranowo, MT.

Insititut Teknologi Bandung 2014

Daftar Isi Penentuan dan Relokasi Hiposenter 1. Grid Search 2. Random Search 3. Simulated Annealing (Guided Random Search) 4. Geiger Jejak Sinar Gelombang (Ray Tracing) & Perhitungan Waktu Tempuh 1. Straight 2. Shooting 3. Bending Gaya Berat Analisis Sinyal 1. Fourier Transform and Inverse Fourier Transform 2. Tapering Function 3. Filter: Highpass, Lowpass, Bandpass 4. Hilbert Transform 5. Autocorrelation, Crosscorrelation, Convolution 6. STFT

Penentuan dan Relokasi Hiposenter Indonesia merupakan negara yang berada pada daerah teknonik aktif. Gempa sudah seharusnya menjadi hal yang biasa terjadi pada daerah ini. Selain gempa karena aktifitas tektonik, gempa juga dapat terjadi dalam skala lokal seperti pada gunung api, gempa mikro karena aktifitas injeksi pada lapangan migas, geothermal, CBM, dan lainnya. Salah satu pengolahan data-data yang telah direkam dalam pengamatan gempa adalah penentuan hiposenter. Hiposenter merupakan posisi sumber gempa dalam 3 dimensi (x,y,z). Banyak metode dan algoritma yang telah berkembang dalam penentuan dan relokasi hiposenter. Dalam pelatihan MATLAB ini kita akan membahas beberapa metode dan algoritma dalam penentuan dan relokasi hiposenter, sebagai berikut : 1. Grid Search 2. Random Search 3. Simulated Annealing (Guided Random Search) 4. Geiger

Grid Search Langkah Membuat data sintetik Diketahui : 1) Terdapat 5 stasiun dan 3 event gempa Koordinat Stasiun : St1 St2 X 8.1212 6.1011 Y 3.7558 1.6615 Z 0.9566 0.1472

St3 7.0149 8.3315 0.8699

St4 0.922 8.3864 0.7694

Event2 9.5169 3.527 5.4093 0.6641

Event3 6.4001 1.8786 5.4635 0.7241

Posisi Hiposenter :

X Y Z T0

Event1 6.2062 2.4733 5.4906 0.0147

2) Kecepatan gelombang P (Vp) = 5.1 km/s

St5 4.2489 4.5161 0.4442

Hitung waktu tempuh (travel time) gelombang gempa dari hiposenter ke stasiun dengan menggunakan rumus : ℎ= / dimana jarak dalam 3-D : = ( − ) +( − ) +( − ) sehingga di dapat waktu observasi : = 0+ ℎ clear,clc % by Rexha Verdhora Ry close all %% MODEL SINTETIK %Menentukan parameter model sintetik xh=1300; %koordinat titik hiposenter yh=1550; zh=-530; vavg=5.1; %velocity rata-rata (km/s) x=[700 2500 1500 1850]; %koordinat stasiun pada sumbu x y=[1900 2000 2500 800]; %koordinat stasiun pada sumbu y z=zeros(1,4); %koordinat stasiun pada sumbu z %Menghitung travel time sintetik pada tiap stasiun for i=1:length(x) tsint(i)=sqrt(((xh-x(i))^2)+((yh-y(i))^2)+((zh-z(i))^2))/vavg; end %Membuat t observasi dengan emberikan noise pada data sintetik noise=0.005; tobs=tsint+normrnd(0,noise*tsint,size(tsint)); %% SYSTEMATICAL GRID SEARCH %Membuat lebar grid L=10; x0=(500:L:2500); y0=(500:L:2500); z0=[0:-L:-600]; P=10; %Perhitungan missfit di setiap grid h = waitbar(0,'Please wait...'); for l=1:length(z0) waitbar(l/length(z0)) for j=1:length(x0) for k=1:length(y0) for i=1:length(x) tcal(i)=sqrt(((x0(j)-x(i))^2)+... ((y0(k)-y(i))^2)+((z0(l)-z(i))^2))/vavg; end dt=tobs-tcal; missfit(k,j)=sqrt(mean(dt.^2)); if missfit(k,j)

end end end close(h) %Plot lokasi hiposenter figure plot3(HypoX,HypoY,HypoZ,'ko','MarkerFaceColor','g','MarkerSize',10) hold on plot3(x,y,z,'v','MarkerFaceColor','k'); title('\fontsize{13}Map of Event Distribution') xlabel('West - East (m)') ylabel('North - South (m)') zlabel('Depth (m)') grid on hold off

Random Search clear,clc % by Rexha Verdhora Ry close all %% MODEL SINTETIK %Menentukan parameter model sintetik xh=1300; %koordinat titik hiposenter yh=1550; zh=-530; vavg=5.1; %velocity rata-rata (km/s) x=[700 2500 1500 1850]; %koordinat stasiun pada sumbu x y=[1900 2000 2500 800]; %koordinat stasiun pada sumbu y z=zeros(1,4); %koordinat stasiun pada sumbu z %Menghitung travel time sintetik pada tiap stasiun for i=1:length(x) tsint(i)=sqrt(((xh-x(i))^2)+((yh-y(i))^2)+((zh-z(i))^2))/vavg; end %Membuat t observasi dengan emberikan noise pada data sintetik noise=0.005; tobs=tsint+normrnd(0,noise*tsint,size(tsint));

%% RANDOM SEARCH %Membuat lebar grid x0=randint(1,5000,[500 2500]); y0=randint(1,5000,[500 2500]); z0=randint(1,5000,[-600 0]); P=10; %Perhitungan missfit di setiap grid h = waitbar(0,'Please wait...'); for l=1:length(x0) waitbar(l/length(x0)) for i=1:length(x) tcal(i)=sqrt(((x0(l)-x(i))^2)+((y0(l)-y(i))^2)... +((z0(l)-z(i))^2))/vavg; end dt=tobs-tcal; missfit(l)=sqrt(mean(dt.^2)); if missfit(l)

Simulated Annealing Flow Chart:

(Grandis,2009)

clear,clc close all

% by Rexha Verdhora Ry

%% MODEL SINTETIK %Menentukan parameter model sintetik xh=1300; %koordinat titik hiposenter yh=2050; zh=-1000; vavg=5.1; %velocity rata-rata (km/s) x=[700 1900 1200 1350]; %koordinat stasiun pada sumbu x y=[2300 2000 3500 800]; %koordinat stasiun pada sumbu y z=zeros(1,4); %koordinat stasiun pada sumbu z %Menghitung travel time sintetik pada tiap stasiun for i=1:length(x) tsint(i)=sqrt(((xh-x(i))^2)+((yh-y(i))^2)+((zh-z(i))^2))/vavg; end %Membuat t observasi dengan emberikan noise pada data sintetik noise=0.005; tobs=tsint+normrnd(0,noise*tsint,size(tsint)); %% INVERSI NON LINIER %Simulated Annealing Inversion Xmin=0; Xmax=3000; Ymin=0; Ymax=3000; Zmin=-1500; Zmax=0;

X=Xmin+rand*(Xmax-Xmin); Y=Ymin+rand*(Ymax-Ymin); Z=Zmin+rand*(Zmax-Zmin); for i=1:length(x) tcal(i)=sqrt(((X-x(i))^2)+((Y-y(i))^2)+((Z-z(i))^2))/vavg; end dt=tobs-tcal; E=dt*dt'; T=E; k=0; U=0; Umax=100; Lmax=100; Amin=100; h = waitbar(0,'Please wait...'); while U<=Umax Lk=0; Ak=0; waitbar(U/Umax) while Lk

if Ak==0 U=U+1; else U=0; end T=T*0.99; k=k+1; end close(h) %Plot kontur dari missfit figure plot3(HypoX,HypoY,HypoZ,'ko',... 'MarkerFaceColor','g','MarkerSize',10) hold on plot3(x,y,z,'v','MarkerFaceColor','k'); title('\fontsize{13}Map of Event Distribution') xlabel('West - East (m)') ylabel('North - South (m)') zlabel('Depth (m)') grid on hold off

Geiger Langkah Metode Geiger 1) Tentukan inisial model (M0) atau posisi hiposenter awal (x0,y0,z0) dan T0. 2) Hitung waktu kalulasi (Tcalci) dari hiposenter tersebut ke stasiun : ( − ) +( − ) +( − = 0+

)

3) Bentuk matriks jacobi (J) dengan ukuran m x 4, dimana m = jumlah stasiun penerima gempa tersebut. ⎡ ⎢ =⎢ ⋮ ⎢ ⎣

⋮

⎤ ⎥ ⋮ ⎥ ⎥ ⎦

⋮

dimana : =

=

=

( (

−

−

) +( (

(

(

−

−

) −

) +(

−

)

) +(

− )

) +(

− )

)

−

) +(

−

(

− )

) +(

−

= 1 4) Hitung waktu deviasi (ΔTi) = Tobsi – Tcalci 5) Lakukan inversi ∆ =( × ) ∆ 6) Kemudian update model (x0,y0,z0,T0) = +∆ 7) Looping tahapan 2) sampai 6) hingga mencapai nilai RMS error sangat kecil =

∑ ΔTi

dimana n = jumlah stasiun yang merekam gempa pada setiap event 8) Lakukan tahapan 1) sampai 7) pada hiposenter gempa berikutnya. 9) Visualisasi hiposenter sintetik. clear,clc close all


%% Bagian 1 %MODEL SINTETIK %Menentukan parameter model sintetik xh=1300; %koordinat titik hiposenter yh=2050; zh=-1000;

vavg=5.1; %velocity rata-rata (km/s) x=[700 1900 1200 1350]; %koordinat stasiun pada sumbu x y=[2300 2000 3500 800]; %koordinat stasiun pada sumbu y z=zeros(1,4); %koordinat stasiun pada sumbu z %Menghitung travel time sintetik pada tiap stasiun for i=1:length(x) tsint(i)=sqrt(((xh-x(i))^2)+((yh-y(i))^2)+((zh-z(i))^2))/vavg; end %Membuat t observasi dengan emberikan noise pada data sintetik noise=0.005; tobs=tsint+normrnd(0,noise*tsint,size(tsint));

%% Bagian 2 %INVERSI NON LINIER %Menentukan parameter model awal x0=700; y0=1500; z0=-700; %Iterasi untuk forward modeling sebanyak 10 kali for j=1:10 %Menghitung travel time kalkulasi dan menyusun matriks Jacobi for i=1:length(x) tcal(i)=sqrt(((x0-x(i))^2)+((y0-y(i))^2)+((z0-z(i))^2))/vavg; dtdx0(i)=(x0-x(i))/(vavg*sqrt(((x0-x(i))^2)+((y0-y(i))^2)+((z0z(i))^2))); dtdy0(i)=(y0-y(i))/(vavg*sqrt(((x0-x(i))^2)+((y0-y(i))^2)+((z0z(i))^2))); dtdz0(i)=(z0-z(i))/(vavg*sqrt(((x0-x(i))^2)+((y0-y(i))^2)+((z0z(i))^2))); end dt=tobs'-tcal'; error(j)=sqrt(mean(dt.^2)); J=[dtdx0' dtdy0' dtdz0']; %Inversi non-linier dm=inv(J'*J)*J'*dt;


figure(2) plot3(x0,y0,z0,'*') hold on grid on %Parameter model baru hasil inversi x0=x0+dm(1); y0=y0+dm(2); z0=z0+dm(3); end HypoX=x0; HypoY=y0; HypoZ=z0;

%Plot posisi stasiun dan hiposenter hasil inversi figure(2) plot3(HypoX,HypoY,HypoZ,'ko',... 'MarkerFaceColor','g','MarkerSize',10) hold on plot3(x,y,z,'v','MarkerFaceColor','k'); title('\fontsize{13}Map of Event Distribution') xlabel('West - East (m)') ylabel('North - South (m)') zlabel('Depth (m)') grid on hold off %Plot error terhadap iterasi figure(1) plot(error) title('RMS Tiap Iterasi') xlabel('Iterasi ke-') ylabel('RMS (s)')

Jejak Sinar Gelombang (Ray Tracing) & Perhitungan Waktu Tempuh Rekonstruksi penjejakan sinar gelombang atau dikenal dengan ray tracing dibutuhkan dalam menghitung waktu tempuh gelombang dari sumber ke stasiun penerima. Ray tracing merupakan bagian penting dalam beberapa metode atau algoritma dalam geofisika, seperti tomografi dan penentuan hiposenter. Dalam modul ini akan dibahas beberapa metode ray tracing, seperti : 1. Straight 2. Shooting (Snell’s Law) 3. Bending (Fermat’s Principle) Kita akan mencoba membangun algoritma sederhana untuk mengaplikasi metodemetode ray tracing di atas. Model kecepatan yang digunakan 2 lapis dimana batas lapisan tegak dengan kecepatan pada lapisan ini 1000 m/s dan 2000 m/s serta posisi source (0,8) dan receiver (10,0). Straight %% Input Data & Parameter % by Ahmad Syahputra Source.X = 0; Source.Y = 8; Receiver.X = 10; Receiver.Y = 0; x0 = 0; y0 = 0; dx = 5; dy = 10; Xgrid = x0:dx:10; Ygrid = y0:dy:10; nx = length(Xgrid)-1; ny = length(Ygrid)-1; [Mesh.X Mesh.Y] =meshgrid(Xgrid,Ygrid); Velocity = [1000 2000]; VelocityPlot = [Velocity; Velocity(size(Velocity,1),:)]; VelocityPlot = [VelocityPlot VelocityPlot(:,size(VelocityPlot,2))]; %% Straight Ray Tracing RayPath.X = linspace(Source.X,Receiver.X,100); RayPath.Y = linspace(Source.Y,Receiver.Y,100); %% Waktu Tempuh Gelombang for i=1:length(RayPath.X) iblok = fix((RayPath.X(i)-x0)/dx)+1; jblok = fix((RayPath.Y(i)-y0)/dy)+1;

nblokpath(i) = (jblok-1)*nx+iblok; if i
Shooting %% Input Data & Parameter % by Ahmad Syahputra Source.X = 0; Source.Y = 8; Receiver.X = 10; Receiver.Y = 0; x0 = 0;y0 = 0; dx = 5;dy = 10; Xgrid = x0:dx:10; Ygrid = y0:dy:10; nx = length(Xgrid)-1; ny = length(Ygrid)-1; [Mesh.X Mesh.Y] =meshgrid(Xgrid,Ygrid); Velocity = [1000 2000]; VelocityPlot = [Velocity; Velocity(size(Velocity,1),:)]; VelocityPlot = [VelocityPlot VelocityPlot(:,size(VelocityPlot,2))]; %% Ray Tracing Shooting SudutDatang = 0:0.1:90; BatasLapisanX = 5; ErrorMin = inf; for i=1:length(SudutDatang) dyy = tan(deg2rad(SudutDatang(i)))*dx; if Source.Y > Receiver.Y BatasLapisanY = Source.Y - dyy; else BatasLapisanY = Source.Y + dyy; end SudutTransmisi = sin(deg2rad(SudutDatang(i)))*Velocity(2)/Velocity(1); dyy = tan(SudutTransmisi)*dx; if Source.Y > Receiver.Y ReceiverTempY = Source.Y - dyy; else ReceiverTempY = Source.Y + dyy; end ErrorDist = sqrt((Receiver.X-Receiver.X)^2+(ReceiverTempYReceiver.Y)^2); if ErrorDist < ErrorMin ErrorMin = ErrorDist; RayPath.X = [Source.X BatasLapisanX Receiver.X]; RayPath.Y = [Source.Y BatasLapisanY ReceiverTempY]; end end %% Time Calculation TimeCalculation = 0; for i=1:length(RayPath.X)-1 TimeCalculation = TimeCalculation + sqrt((RayPath.X(i)RayPath.X(i+1))^2+(RayPath.Y(i)-RayPath.Y(i+1))^2)/Velocity(i); end %% Visualisasi pcolor(Mesh.X,Mesh.Y,VelocityPlot)

set(gca,'Ydir','Reverse') hold on plot(Source.X,Source.Y,'*r') hold on plot(Receiver.X,Receiver.Y,'vb') hold on plot(RayPath.X,RayPath.Y,'c') hold off title('Shooting Ray Tracing') xlabel('Distance') ylabel('Depth') colorbar

Bending %% Input Data & Parameter % by Ahmad Syahputra Source.X = 0; Source.Y = 8; Receiver.X = 10; Receiver.Y = 0; x0 = 0;y0 = 0; dx = 5;dy = 10; Xgrid = x0:dx:10; Ygrid = y0:dy:10; nx = length(Xgrid)-1; ny = length(Ygrid)-1; [Mesh.X Mesh.Y] =meshgrid(Xgrid,Ygrid); Velocity = [1000 2000]; VelocityPlot = [Velocity; Velocity(size(Velocity,1),:)]; VelocityPlot = [VelocityPlot VelocityPlot(:,size(VelocityPlot,2))];

%% Ray Tracing Bending %Bending Point BatasLapisanY = 0:1:10; BatasLapisanX = 5*ones(1,length(BatasLapisanY)); minTime = inf; for i=1:length(BatasLapisanY) Time1 = sqrt((Source.X-BatasLapisanX(i))^2+(Source.YBatasLapisanY(i))^2)/Velocity(1); Time2 = sqrt((Receiver.X-BatasLapisanX(i))^2+(Receiver.YBatasLapisanY(i))^2)/Velocity(2); TimeCalculation = Time1+Time2; if TimeCalculation < minTime minTime = TimeCalculation; RayPath.X = [Source.X BatasLapisanX(i) Receiver.X]; RayPath.Y = [Source.Y BatasLapisanY(i) Receiver.Y]; end end %% Visualisasi pcolor(Mesh.X,Mesh.Y,VelocityPlot) set(gca,'Ydir','Reverse') hold on plot(Source.X,Source.Y,'*r') hold on plot(Receiver.X,Receiver.Y,'vb') hold on plot(RayPath.X,RayPath.Y,'c') hold on plot(BatasLapisanX,BatasLapisanY,'.w') title('Bending Ray Tracing') xlabel('Distance') ylabel('Depth') colorbar

Gaya Berat Inversi Non-Linier Anomaly Bola clear,clc close

% Rexha Verdhora Ry

%% MODEL SINTETIK %Menentukan parameter model sintetik G=6.67*10^(-11); r=150; noise=0.03; rho=2500; xbola=400; zbola=650; %Menentukan titik pengukuran x=linspace(0,1000,21); z=zeros(1,21); %Menghitung nilai g sintetik di setiap titik pengukuran for i=1:21 a=((x(i)-xbola)^2+(z(i)-zbola)^2)^(3/2); g(i)=G*(4/3)*pi*(r^3)*zbola*rho/a; end %Membuat g sintetik dalam mgal gsint=10^5*g'; %Membuat plot dari data g sintetik figure(1) subplot(2,3,1) plot(x,gsint,'.') title('Data Sintetik') xlabel('Jarak') ylabel('Gaya berat (mgal)') %Membuat plot model sintetik dari bola z1bola=-zbola;t=linspace(0,2*pi,200); xh=r*cos(t)+xbola;zh=r*sin(t)+z1bola; subplot(2,3,4) plot(xh,zh,'-b') title('Model Sintentik') xlabel('Jarak');ylabel('Kedalaman') axis ([0 1000 -1000 0]) grid on %Membuat g observasi dengan menambahkan noise thd model sintetik gobs=gsint+normrnd(0,noise*gsint,size(gsint)); %Membuat plot dari data g observasi subplot(2,3,2) plot(x,gobs,'r.') title('Data Sintetik + noise') xlabel('Jarak') ylabel('Gaya berat (mgal)')

%% INVERSI NON LINIER %Menentukan parameter awal x0=150; z0=150; %Iterasi untuk forward modeling sebanyak 20 kali for j=1:20 %Membuat plot model bola dari forward modeling z01=-z0;t0=linspace(0,2*pi,200); xh0=r*cos(t)+x0;zh0=r*sin(t)+z01; subplot(2,3,5) plot(xh0,zh0,'-') title('Model per Iterasi') xlabel('Jarak');ylabel('Kedalaman') hold on grid on axis ([0 1000 -1000 0]) %Menghitung g forward modeling for i=1:21 a=((x(i)-x0)^2+(z(i)-z0)^2)^(3/2); g1(i)=G*(4/3)*pi*(r^3)*z0*rho/a; end %Membuat g forward modeling dalam mgal gfm=10^5*g1'; d=gobs-gfm;

%selisih g observasi dengan g kalkulasi

%Membuat matriks jacobbi for i=1:21 a=((x(i)-x0)^2+(z(i)-z0)^2)^(5/2); jx0(i)=G*(4/3)*pi*(r^3)*z0*rho*2*(x(i)-x0)/a; jz0(i)=G*(4/3)*pi*(r^3)*z0*rho*2*(z(i)-z0)/a; end J(:,1)=jx0'; J(:,2)=jz0'; %Inversi non linier dm=inv(J'*J)*J'*d; %Membuat matriks dari parameter model awal m=[x0 z0]'; %Rexha Verdhora Ry %Menghitung parameter model baru hasil inversi %10^-5 sebagai konversi ke meter m=m+(dm*10^(-5)); x0=m(1); z0=m(2); missfit(j)=sqrt((gfm-gobs)'*(gfm-gobs)); end

%Membuat plot model bola dari forward modeling subplot(2,3,5) plot(xh0,zh0,'r-') title('Model per Iterasi') xlabel('Jarak') ylabel('Kedalaman') hold on grid on axis ([0 1000 -1000 0]) for i=1:21 a=((x(i)-x0)^2+(z(i)-z0)^2)^(3/2); g1(i)=G*(4/3)*pi*(r^3)*z0*rho/a; end gfm=10^5*g1' %Membuat plot dari fit g observasi dengan g forward modeling subplot(2,3,3) plot(x,gobs,'r.') hold on plot(x,gfm,'-k') title('Data Sintetik + Fit Forward Modeling') xlabel('Jarak') ylabel('Gaya berat (mgal)') %Membuat plot model bola yang paling fit dari forward modeling subplot(2,3,6) plot(xh0,zh0,'r-') title('Model hasil inversi') xlabel('Jarak') ylabel('Kedalaman') grid on axis ([0 1000 -1000 0])

Analisis Sinyal Fourier Transform and Inverse Fourier Transform Fourier Transform adalah transformasi sinyal dari domain waktu ke domain frekuensi atau dari domain jarak ke domain bilangan gelombang. Untuk mengubah suatu sinyal dari domain waktu ke domain frekuensi dapat dilakukan di MATLAB dengan cara: Xf = fft(xt)

Sedangkan untuk membalikannya ke domain frekuensi dilakukan Inverse Fourier Transform: xt = ifft(Xf)

Dimana adalah sinyal pada domain waktu dan adalah sinyal dalam domain frekuensi. Pada Fourier Transform, hasilnya akan terbagi 2 yaitu bilangan real dan imaginer dalam bentuk + dimana

adalah bilangan real,

adalah bilangan imaginer, dan adalah unit imaginer (√−1).

Contoh soal: % by Waskito Pranowo dt = 0.002; % sampling rate in second t = 0:dt:1; % time xt = sin(2*pi*80*t) + cos(2*pi*40*t); % signal in time domain Xf = fft(xt)*dt; % signal in frequency domain X = Xf(1:round(length(Xf)/2)); % single-sided spectrum f = linspace(0,1,length(X))/(dt*2); % frequency Xreal = real(X); % real part of signal in frequency domain Ximag = imag(X); % imaginary part of signal in frequency domain Xamp = abs(X); % amplitude of signal in frequency domain figure('Color','w') subplot(2,2,1) plot(t,xt) xlabel('Time (second)'), title('Signal in time domain') subplot(2,2,2) plot(f,Xamp) xlabel('Frequency (Hz)'), title('Amplitude Spectrum') subplot(2,2,3) plot(f,Xreal) xlabel('Frequency (Hz)'), title('Real Part') subplot(2,2,4) plot(f,Ximag) xlabel('Frequency (Hz)'), title('Imaginary Part')

Tapering Function atau Apodization Function Tapering function adalah sebuah fungsi yang digunakan untuk menjadikan nilai sinyal menurun ke nol secara bertahap pada daerah tepi sinyal.

Beberapa contoh window/ tapering function (Sumber: Weisstein, Eric W. "Apodization Function." From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/ApodizationFunction.html)

Jenis-jenis filter dapat dilihat di Help  search: window.

Lowpass, Bandpass, dan Highpass Filter Filter digunakan untuk mengambil bagian frekuensi yang kita inginkan dan membuang bagian frekuensi yang tidak diinginkan.   

Lowpass filter, meloloskan sinyal berfrekuensi rendah dan membuang sinyal berfrekuensi tinggi. Highpass filter, meloloskan sinyal berfrekuensi tinggi dan membuang sinyal berfrekuensi rendah. Bandpass filter, meloloskan frekuensi dalam rentang tertentu dan membuang frekuensi diluar rentang tersebut.

Script untuk membangun filter menggunakan MATLAB:

clear,clc

% by Waskito Pranowo

dt = 0.001; % sampling rate (in second) fc = [10 20 80 100]; %Cut-off Frequency %Filter Builder f1 = fc(1);f2 = fc(2);f3 = fc(3);f4 = fc(4); y = [zeros(1,f1) linspace(0,1,(f2-f1)+1)... ones(1,f3-f2-1) linspace(1,0,(f4-f3)+1)... zeros(1,1000/(2*dt*1000)-f4-1)]; F = linspace(0,1,length(y))/(2*dt); figure plot(F,y) title('Filter in Frequncy Domain') %Inverse Fast Fourier Transform y2 = [y fliplr(y(2:end))]; n=length(y2); Y = real(ifft(y2)); Y = [Y(round(n/2)+1:end) Y(1:round(n/2))]; t = (-floor(n/2):1:floor(n/2))*dt; figure plot(t,Y) title('Filter in Time Domain')

Lowpass-Filter dengan dengan cut-off frekuensi 0 dan 70 Hz

Highpass-Filter dengan dengan cut-off frekuensi 100 Hz

Bandpass-Filter dengan dengan cut-off frekuensi 10-20-80-100 Hz

Contoh soal: clear,clc


dt = 0.001; % sampling rate in second t = 0:dt:1; % time xt = sin(2*pi*90*t) + cos(2*pi*30*t); % signal in time domain

fc = [0 0 45 70]; %Cut-off Frequency %Filter Builder f1 = fc(1);f2 = fc(2);f3 = fc(3);f4 = fc(4); y = [zeros(1,f1) linspace(0,1,(f2-f1)+1)... ones(1,f3-f2-1) linspace(1,0,(f4-f3)+1)... zeros(1,1000/(2*dt*1000)-f4-1)]; F = linspace(0,1,length(y))/(2*dt); % Frequency figure('Color','w') plot(F,y) title('Filter in Frequncy Domain') y2 = [y fliplr(y(2:end))]; n=length(y2); % double-side Filter Ff = linspace(0,1,length(y2))/(dt); Ff2 = linspace(0,1,length(xt))/(dt); Y= interp1(Ff,y2,Ff2); Xf = fft(xt,length(Y))*dt; % signal in frequency domain X = Xf(1:round(length(Xf)/2)); % single-sided spectrum FiltX = Xf.*Y; % Filter the signal Xt = real(ifft(FiltX))/dt; % Filtered signal in time domain figure('Color','w') subplot(2,1,1) plot(t,xt) subplot(2,1,2) plot(t,Xt)

Hilbert Transform Hilbert transform digunakan untuk merotasi suatu sinyal sejauh 90°. Signal asli disebut sinyal real, sedangkan yang telah dirotasikan 90° disebut sebagai quadrature. y = hilbert(x) q = -imag(y);

Dimana adalah sinyal real, adalah Hilbert transform dan adalah quadrature. Pada Hilbert Transform, hasilnya akan terbagi 2 yaitu bilangan real dan quadrature dalam bentuk + dimana adalah bilangan real, adalah bilangan quadrature, dan adalah unit imaginer (√−1). Walaupun hasil dari Hilbert Transform memiliki format yang mirip dengan Fourier Transform, tetapi dan di sini memiliki pengertian yang berbeda dengan Fourier Transform. clear, clc


dt = 0.001; t = 0:dt:0.5; x y q A

= = = =

cos(2*pi*20*t); % Real signal hilbert(x); % Hilbert transform -imag(y); % Quadrature signal (real signal yang dirotasi 90 derajat) abs(y); % Instantaneous Amplitude

figure('Color','w') subplot(3,1,1) plot(t,x) xlabel('Time (s)'), ylabel('Amplitude'), title('Real Signal') ylim([-1.1 1.1]) subplot(3,1,2) plot(t,q) xlabel('Time (s)'), ylabel('Amplitude'), title('Quadrature Signal') ylim([-1.1 1.1]) subplot(3,1,3) plot(t,A) xlabel('Time (s)'), ylabel('Amplitude'), title('Instantaneous Amplitude')

Convolution and Correlation Konvolusi adalah proses penggabungan dua sinyal dimana prosesnya meliputi: 1. 2. 3. 4. 5.

Salah satu sinyal diputar atau dicerminkan, Kedua sinyal dikalikan, Hasil perkalian dijumlahkan, Sinyal yang telah diputar pada nomor 1 digeser, Kembali ke nomor 2.

Sedangkan proses korelasi mirip dengan konvolusi tetapi tidak mengalami tahap pemutaran (tahap nomor 1 dalam konvolusi). clear,clc x_1 = [1 3 5 2 7]; x_2 = [1 4 2 4 5]; % Konvolusi Y1 = conv(x_1,x_2) Y2 = conv(x_2,x_1) % Korelasi Z1 = xcorr(x_1,x_2) Z2 = xcorr(x_2,x_1)

Hasil dari proses tersebut adalah: Y1 = 1

7 19 32 42 67 47 38 35

Y2 = 1

7 19 32 42 67 47 38 35

Z1 = 5.0000 19.0000 39.0000 40.0000 66.0000 55.0000 27.0000 30.0000 7.0000 Z2 = 7.0000 30.0000 27.0000 55.0000 66.0000 40.0000 39.0000 19.0000 5.0000

Dapat dilihat bahwa pada proses konvolusi berlaku sifat komutatif, proses korelasi tidak berlaku sifat komutatif,  ≠  .

∗

=

∗ , sedangkan pada

Short Time Fourier Transform (STFT) 1946, Dennis Gabor mengadaptasi Fourier transform untuk menganalisis bagian kecil suatu sinyal pada waktu tertentu. Teknik ini disebut windowing sinyal. Adaptasi Gabor, disebut ShortTime Fourier Transform (STFT), memetakan sinyal menjadi 2D antara waktu dan frekuensi. (Matlab Help, 2008)

Konsep STFT (Matlab Help, 2008)

Pada STFT, jendela yang digunakan berbentuk boxcar.

∞

( ) ( − )

( , )= ∞

dimana

( ) adalah fungsi pada domain waktu dan

( ) adalah window. Window yang

digunakan dapat berbentuk boxcar atau tapering function lainnya. Gabor transform adalah ketika STFT menggunakan Gaussian Window. Aplikasi STFT dan Gabor Transform pada tugas ini diterapkan pada fungsi berikut, cos 2 cos 2 ( )= cos 2 cos 2 clear,clc

10 , 50 , 100 , 150 ,

0.0 ≤ < 0.5 0.5 ≤ < 1.0 1.0 ≤ < 1.5 1.5 ≤ < 2.0


dt = 0.002; t = 0:dt:2; for i = 1:length(t) if t(i) < 0.5 x(i) = cos(2*pi*10*t(i)); elseif t(i) >= 0.5 & t(i) < 1 x(i) = cos(2*pi*50*t(i)); elseif t(i) >= 1 & t(i) < 1.5 x(i) = cos(2*pi*100*t(i)); else x(i) = cos(2*pi*150*t(i)); end end figure('Color','w') plot(t,x) xlabel('Time (s)'), ylabel('Amplitude'), title('Signal in Time Domain') [S,F,T] = spectrogram(x,100,90,length(x),1/dt); figure('Color','w') imagesc(T,F,abs(S)) xlabel('Time (s)'), ylabel('Frequency (Hz)')

Modul Pelatihan Matlab 2014

Recommend Documents