fudaopu.m

来自「关于功率谱估计的,周期图法,自相关法,bartlett法和welch法,还有一个」· M 代码 · 共 53 行

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('最小相位法的复倒谱');