function phwave()
# read file
 load exp.d
 t=exp(:,1);
 d=exp(:,2);

#draw wave form
 n=16
 dt=0.5
 t=dt*(0:n-1);
 subplot(2,1,1);
 plot(t,d);
 grid on;

#fft
 tl=dt*n
 df=1/tl
 f=df*(0:n-1);
 fftd=fft(d)
 afd=abs(fftd);

#inverse fft
 fftpd=fftd./afd;
 fftpd(1,1)=0;
 fftpd
 pd=ifft(fftpd);
 pd
 subplot(2,1,2);
 plot(t,pd);
 grid on;

endfunction