1 同態處理

訊號的同態處理也稱同態濾波。大概步驟為：

f(x,y)→ln→DFT→H(u,v)→(DFT)-1→exp→g(x,y)

雖然，一般用於影象處理。但是，博主將同態濾波用於語音訊號的濾波。直接上程式吧

clc;clear
%%
filedir=[];                             % 指定檔案路徑
filename='bluesky3.wav';                % 指定檔名
fle=[filedir filename]                  % 構成路徑和檔名的字串
[xx,fs]=audioread(fle);                   % 讀入資料檔案
xx=xx-mean(xx);                         % 消除直流分量
x=xx/max(abs(xx));                      % 幅值歸一化
SNR=5;                                  % 設定SNR
I=Gnoisegen(x,SNR);                % 疊加噪聲
N=length(x);                            % 訊號長度
time=(0:N-1)/fs;                        % 設定時間

%%
%同態濾波
[M,N]=size(I);  
rL=0.2;  
rH=7.7;%可根據需要效果調整引數  
c=2;  
d0=10;  
I1=log(I+1);%取對數  
FI=fft2(I1);%傅立葉變換  
n1=floor(M/2);  
n2=floor(N/2);  
for i=1:M  
    for j=1:N  
        D(i,j)=((i-n1).^2+(j-n2).^2);  
        H(i,j)=(rH-rL).*(exp(c*(-D(i,j)./(d0^2))))+rL;%高斯同態濾波  
    end  
end  
I2=ifft2(H.*FI);%傅立葉逆變換  
output=real(exp(I2))-1; 
%%
% 作圖
subplot 311; plot(time,xx,'k'); ylabel('幅值') 
ylim([-1 1 ]); title('原始語音訊號');
subplot 312; plot(time,I,'k'); ylabel('幅值') 
ylim([-1 1 ]); title('帶噪語音訊號');
subplot 313; plot(time,output,'k'); 
ylim([-1 1 ]); title('同態濾波輸出語音訊號');
xlabel('時間/s'); ylabel('幅值')
 

結果

       當語音訊號加入信噪比為5的高斯白噪聲之後，得到了帶噪語音訊號，經過同態濾波之後，雖然達不到完全濾波，但是跟帶噪語音相比，達到了降噪效果。這裡，博主並未對語音訊號進行分幀，之前試過，效果都差不多，最後還需要把分幀後的訊號進行整合。語音訊號與影象訊號不同之處在於前者是一維，後者是二維。博主之前從電影上面剪下一段語音，來進行語音處理，但是發現無法處理。這是因為剪下來的語音的雙聲道，而我們一般處理的語音訊號都是單聲道，所以博主最後用Goldwave軟體，將雙聲道語音變為單聲道。

2 倒譜分析

                                                                                            （1）

序列c(n)是x(n)對數幅值譜的逆傅立葉變換，稱做倒頻譜，簡稱倒譜。這裡，c(n)的量綱是Quefrency，又被稱做倒頻，它實際的單位還是時間的單位。

       例子：從su1.txt中讀入語音訊號，訊號取樣頻率是16000Hz，要求如下：

a.按照公式（1）計算訊號的倒譜；

b.把語音的聲門激勵訊號和聲道衝激響應分離，分別計算相應的頻譜。

clear all; clc; close all;
y=load('su1.txt');                            % 讀入資料
fs=16000; nfft=1024;                          % 取樣頻率和FFT的長度
time=(0:nfft-1)/fs;                           % 時間刻度
figure(1), subplot 211; plot(time,y,'k');     % 畫出訊號波形
title('訊號波形'); axis([0 max(time) -0.7 0.7]);
ylabel('幅值'); xlabel(['時間/s' 10 '(a)']); grid;
figure(2)
nn=1:nfft/2; ff=(nn-1)*fs/nfft;               % 計算頻率刻度
Y=log(abs(fft(y)));                           % 按式(1)取實數部分
subplot 211; plot(ff,Y(nn),'k'); hold on;     % 畫出訊號的頻譜圖
z=ifft(Y);                                    % 按式(1)求取倒譜
figure(1), subplot 212; plot(time,z,'k');     % 畫出倒譜圖
title('訊號倒譜圖'); axis([0 time(512) -0.2 0.2]); grid; 
ylabel('幅值'); xlabel(['倒頻率/s' 10 '(b)']);
mcep=29;                                      % 分離聲門激勵脈衝和聲道衝激響應
zy=z(1:mcep+1);
zy=[zy' zeros(1,1000-2*mcep-1) zy(end:-1:2)']; % 構建聲道衝激響應的倒譜序列
ZY=fft(zy);                                   % 計算聲道衝激響應的頻譜
figure(2),                                    % 畫出聲道衝激響應的頻譜，用灰線表示
line(ff,real(ZY(nn)),'color',[.6 .6 .6],'linewidth',3);
grid; hold off; ylim([-4 5]);
title('訊號頻譜（黑線）和聲道衝激響頻譜（灰線）')
ylabel('幅值'); xlabel(['頻率/Hz' 10 '(a)']); 

ft=[zeros(1,mcep+1) z(mcep+2:end-mcep)' zeros(1,mcep)]; % 構建聲門激勵脈衝的倒譜序列
FT=fft(ft);                                  % 計算聲門激勵脈衝的頻譜
subplot 212; plot(ff,real(FT(nn)),'k'); grid;% 畫出聲門激勵脈衝的頻譜
title('聲門激勵脈衝頻譜')
ylabel('幅值'); xlabel(['頻率/Hz' 10 '(b)']); 
 

結果如下：

這裡，博主想重點說一下以下三行程式。

mcep=29;

zy=[zy' zeros(1,1000-2*mcep-1) zy(end:-1:2)']; % 構建聲道衝激響應的倒譜序列（預加重）

ft=[zeros(1,mcep+1) z(mcep+2:end-mcep)' zeros(1,mcep)]; % 構建聲門激勵脈衝的倒譜序列（預加重）

 在這個例子中，取樣頻率為16000HZ，認為基音訊率都應低於500HZ，為了留有餘量，取16000/550=29，即在倒譜域中第29條譜線之前是反映了包絡的係數。所以取29就是這樣得來的。   訊號模型可看成激勵模型、聲道模型、輻射模型的串聯。預加重目的是為了更好的分析聲道，去掉激勵模型和輻射模型的影響。激勵模型是一個二階的低通模型，輻射模型是一階高通模型，所以再增加一個預加重，增加加一個一階高通，以平衡激勵模型和輻射模型的影響。   倒譜中橫座標不是頻率，而是倒頻率（quefrency，量綱是時間），它的刻度同時間序列中x對應的刻度相同，在FFT(x)以後包絡的起伏比較緩，就是它的倒頻率低。倒頻譜中的橫軸也反映了在取樣頻率16000下不同頻率的週期，550Hz的週期為29，即第29條譜線。在大於第29條譜線以後的譜線，對應的訊號分量的頻率都要小於550Hz，即在基頻的區間中。在倒頻率範 圍中我們取小於550Hz的倒頻譜來求包絡。  3 Mel頻率倒譜系數的分析 這裡，博主覺得理解了兩個知識點，Mel係數就算弄明白了。 3.1 知識點1 Mel濾波器組 程式呼叫的函式如下： x=melbankm(24,256,8000,0,0.5,w) 這裡，訊號取樣頻率為8000Hz，每幀長256，設定為24個mel濾波器組。 涉及的知識點
這裡，f(m)為中心頻率。fl為濾波器頻率範圍內的最低頻率；fh為最高；N為DFT（FFT）時的長度；fs為取樣頻率。最後，每個帶通濾波器的傳遞函式為
知識點2： 對語音訊號進行如下操作，最終可以得到Mel濾波器的能量S(i,m) 預加重→分幀→加窗→倒譜
把Mel濾波器的能量取對數後計算DCT：
這裡，S是Mel濾波器能量；m是第m個Mel濾波器（共有M個）；i是第i幀；n是DCT的譜線。

function ccc=mfcc_m(x,fs,p,frameSize,inc)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                 function ccc=mfcc_m(x);
%對輸入的語音序列x進行MFCC引數的提取，返回MFCC引數和一階
%差分MFCC引數，Mel濾波器的個數為p，取樣頻率為fs
%對x每frameSize點分為一幀，相鄰兩幀之間的幀移為inc
% fft變換的長度為幀長
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 按幀長為frameSize，Mel濾波器的個數為p，取樣頻率為fs
% 提取Mel濾波器引數，用海明窗函式
bank=melbankm(p,frameSize,fs,0,0.5,'m');
% 歸一化mel濾波器組係數
bank=full(bank);
bank=bank/max(bank(:));

% DCT係數,12*p
for k=1:12
  n=0:p-1;
  dctcoef(k,:)=cos((2*n+1)*k*pi/(2*p));
end

% 歸一化倒譜提升視窗
w = 1 + 6 * sin(pi * [1:12] ./ 12);
w = w/max(w);

% 預加重濾波器
xx=double(x);
xx=filter([1 -0.9375],1,xx);

% 語音訊號分幀
xx=enframe(xx,frameSize,inc);
n2=fix(frameSize/2)+1;
% 計算每幀的MFCC引數
for i=1:size(xx,1)
  y = xx(i,:);
  s = y' .* hamming(frameSize);
  t = abs(fft(s));
  t = t.^2;
  c1=dctcoef * log(bank * t(1:n2));
  c2 = c1.*w';
  m(i,:)=c2';
end

%差分系數
dtm = zeros(size(m));
for i=3:size(m,1)-2
  dtm(i,:) = -2*m(i-2,:) - m(i-1,:) + m(i+1,:) + 2*m(i+2,:);
end
dtm = dtm / 3;
%合併mfcc引數和一階差分mfcc引數
ccc = [m dtm];
%去除首尾兩幀，因為這兩幀的一階差分引數為0
ccc = ccc(3:size(m,1)-2,:);

例子：這裡有兩個語音，一個是母音i，另外一個是a。利用MFCC引數觀察它們的匹配情況，計算它們的MFCC距離。這裡，距離我們用馬氏距離來進行計算。

%
% pr3_3_2 
clear all; clc; close all;

[x1,fs]=audioread('s1.wav');      % 讀入訊號s1-\i1\
x3=audioread('a1.wav');           % 讀入訊號a1-\a1\
wlen=200;                       % 幀長
inc=80;                         % 幀移
x1=x1/max(abs(x1));             % 幅值歸一化
x3=x3/max(abs(x3));

% 計算/i1/與/a1/之間的匹配比較
[Dcep,Ccep1,Ccep2]=mel_dist(x1,x3,fs,16,wlen,inc);
figure(2)
plot(Ccep1(3,:),Ccep2(3,:),'k+'); hold on
plot(Ccep1(7,:),Ccep2(7,:),'kx'); 
plot(Ccep1(12,:),Ccep2(12,:),'k^');
plot(Ccep1(16,:),Ccep2(16,:),'kh'); 
legend('第3幀','第7幀','第12幀','第16幀')
xlabel('訊號x1');ylabel('訊號x3')
axis([-12 12 -12 12]);
line([-12 12],[-12 12],'color','k','linestyle','--');
title('/i1/與/a1/之間的MFCC引數匹配比較')

結果如下：

PS:樓主這裡matlab版本是2016b，這個版本里面很多函式與之前的不一樣。比如：wavread變成了audioread。