clear all
%輸入x序列
x=input('請輸入x序列：');
N=length(x);
n=2^(ceil(log2(N)));

%求復倒譜
X=fft(x,n);%對x序列進行離散傅立葉變換，用FFT方式
Xn=log(X+eps);%取自然對數，＋eps避免log0出現
%對x序列進行逆FFT變換
xn=ifft(Xn)
disp('x的復倒譜為：');
disp(xn);

%最小相位法求復倒譜
for i=0:(n-1)
   if (i>=1 && i<=N/2)
     f=2;
   end
   if (i==0 || i==N/2)
     f=1;
   end
   if (i<=N-1 && i>N/2)   
     f=0;
   end
   t1(i+1)=xn(i+1)*f;   
end

%由xn對原x序列進行還原
yn=xn;%對yn序列進行FFT變換
Yn=fft(yn);
Y=exp(Yn);%做指數運算
y=ifft(Y);
disp('還原序列為');
disp(y);

subplot(3,1,1);
stem(x);
title('原序列');

subplot(3,1,2);
stem(y);
title('還原序列');

subplot(3,1,3);
stem(t1);
title('最小相位法的復倒譜');