[信号処理] gwt cepstrum

gwt_cepstrum.py

function [ Ceps_env, Ceps_det ] = Analys_Cepstrum( arg_GWT, arg_lifter )
%   アルゴリズム
% http://jp.mathworks.com/help/signal/ref/rceps.html#bqrvf7e-1
% rceps は、参考文献 [2] のアルゴリズム 7.2 を実装したもので、以下の操作を行います。
% y = real(ifft(log(abs(fft(x)))));
% ケプストラム領域で適切なウィンドウ処理を適用すると、復元された最小位相信号が作成されます。
%
% w = [1;2*ones(n/2-1,1);ones(1-rem(n,2),1);zeros(n/2-1,1)];
% ym = real(ifft(exp(fft(w.*y))));
%
%   h = fft(x);
% logh = log(abs(h)) + sqrt(-1)*rcunwrap(angle(h));
% y = real(ifft(logh));
%
% http://media.sys.wakayama-u.ac.jp/kawahara-lab/LOCAL/diss/diss7/S3_6.htm#6-3
% %------野原アルゴリズム-------
% lifter=3
% %---------FFT--------------------------
% data_fft = abs(fft(data0));
% %---------Cepstrum--------------------------
% cps = real(ifft(log(abs(fft(data0)))));
% %---------Liftering (LowPass)--------------------------
% cps_envelope = cps; %Copy
% % cps_envelope(lifter+1:length(cps)-lifter)=0; 
% lif = [ones(1,lifter), zeros(1, length(cps_envelope)-2*lifter), ones(1,lifter)]';
% cps_envelope = cps_envelope.*lif;
% %---------Liftering (HighPass)--------------------------
% cps_detail = cps; %Copy
% % cps_detail(1 : lifter) = 0;
% % cps_detail(length(cps_detail) - (lifter)+1 : end) = 0;
% lif = [zeros(1,lifter), ones(1, length(cps_detail)-2*lifter), zeros(1,lifter)]';
% cps_detail = cps_detail.*lif;
% %---------iCepstrum--------------------------
% spectrum_envelope = abs(exp(fft(cps_envelope)));
% spectrum_detail = abs(exp(fft(cps_detail)));
% %---------検算--------------------------
% fft_minus_cps = exp(log(data_fft) - log(spectrum_envelope)); %linear
% error = sum(fft_minus_cps - spectrum_detail)
% %---------Figure plot--------------------------
% figure;
% hold on;
% plot(20*log10(data_fft));
% plot(20*log10(spectrum_envelope));
% plot(20*log10(spectrum_detail));
% plot(20*log10(fft_minus_cps));
% hold off;
% legend('fft','env','det','minus');

% arg_GWT: 振幅スペクトル

% リフタリング次数
% lifter_lowpass = arg_lifter;
% lifter_highpass = arg_lifter;

% データの呼び出し
gwtl = arg_GWT;

% 折り返し成分
gwtr = fliplr(gwtl);
gwt = [gwtl, gwtr];

% バッファの格納
FFT_abs=zeros(size(gwt));
Ceps_ceps = zeros(size(gwt));
Ceps_ceps_env = zeros(size(gwt));
Ceps_ceps_det = zeros(size(gwt));
Ceps_sspec=zeros(size(gwt));
Ceps_sspec_env=zeros(size(gwt));
Ceps_sspec_det=zeros(size(gwt));

% 周波数領域
FFT_abs = gwt;

% ケプストラム
% 'Cepstrum'
for count1 = 1:size(FFT_abs,1)
  absSpec = FFT_abs(count1,:);
  Ceps_ceps(count1,:) = real(ifft(log(absSpec)));
end

% CMN
%   Ceps_ceps = detrend(Ceps_ceps,'constant');

% リフタリング
% lifter1 = lifter_lowpass;
% lifter2 = lifter_highpass;
lifter = arg_lifter;
lifter1 = arg_lifter;
lifter2 = arg_lifter;

% バッファのコピー
Ceps_ceps_env = Ceps_ceps;
Ceps_ceps_det = Ceps_ceps;

% ローパス
Ceps_ceps_env( :,  lifter1 : size(Ceps_ceps,2)-(lifter1-1) ) = 0; %間違い
% Ceps_ceps_env( :,  lifter1+1 : size(Ceps_ceps,2)-lifter1 ) = 0; % 修正 野原2015/10/07
% lift_low = [ones(1,lifter), zeros(1, length(Ceps_ceps)-2*lifter), ones(1,lifter)];
% for i = 1:size(Ceps_ceps,1)
%   Ceps_ceps_env(i,:) = Ceps_ceps .* lift_low;
% end

% ハイパス
Ceps_ceps_det( :,  1 : lifter2) = 0;
Ceps_ceps_det( :,  size(Ceps_ceps,2) - (lifter2)+1 :end) = 0;
% lif_high= [zeros(1,lifter), ones(1, length(Ceps_ceps)-2*lifter), zeros(1,lifter)];
% for i = 1:size(Ceps_ceps,1)
%   Ceps_ceps_det(i,:) = Ceps_ceps .* lift_low;
% end

% 逆ケプストラム
% 'iCepstrum'
for count2 = 1:size(FFT_abs,1)
%   Ceps_sspec(count2,:) = abs(exp(fft(Ceps_ceps(count2,:))));
  Ceps_sspec_env(count2,:) = abs(exp(fft(Ceps_ceps_env(count2,:))));
  Ceps_sspec_det(count2,:) = abs(exp(fft(Ceps_ceps_det(count2,:))));
end


% 反転成分の除去
% FFT_abs(:, size(FFT_abs,2)/2+1:end) = [];
% Ceps_sspec(:, size(Ceps_sspec,2)/2+1:end) = [];
Ceps_sspec_env(:, size(Ceps_sspec_env,2)/2+1:end) = [];
Ceps_sspec_det(:, size(Ceps_sspec_det,2)/2+1:end) = [];

Ceps_env = Ceps_sspec_env;
Ceps_det = Ceps_sspec_det;

end