文章主要介紹噪聲消除，文章來自博客園RTC.Blacker，支持原創，轉載必須說明出處，歡迎關注微信公眾號blacker，更多詳見www.rtc.help
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
RTC中聲音處理是個很麻煩的事，難點很多，回聲，噪聲，嘯聲，增益，等等，其實從中國潛艇噪聲那么大就可看出這東西確實不好處理。
這些年我也被這些問題搞得煩死了，特別是android上面的，當然折騰過程中也發現了很多被大家誤解的現象：
1、只有qq和微信語音效果才處理得比較好？
不是，有些公司語音效果也很牛b，只不過他們的產品的終端用戶不是個人，所以我們很少聽說，但這不影響他們在業界的口碑，為了避嫌，名字我就不說了，大家以后有用到的時候自然會知道。
2、做圖像識別和語音識別的公司在RTC方面就一定很牛？
不一定，昨天科技園一家做通訊的公司給我看了他們的一個合作客戶——硅谷很牛的一家做機器人公司（語音識別和視頻交互，不是那種簡單的攝像頭轉一下，掃掃地的），人家對RTC基本上就是一片空白，這也正說明了術業有專攻。
3、視頻通訊的都是用WebRTC搞的？
有些公司是完全用WebRTC搞的（用libjingle，ICE，webrtc），而更多公司只是用他里面的一些模塊，當然也有公司基本上都是自己搞的（包括編解碼，傳輸，抖動緩沖，回聲消除，噪聲消除等），WebRTC對他們來說主要是做參考。
4、找個大牛用webrtc就可以搞定RTC所有問題？
用之前在微信里面做底層編解碼的朋友的話解答：“我們做底層算法的不熟悉上層應用，熟悉上層應用的不熟悉底層算法”，我當時補充了一句：“兩者都能搞定的不熟悉市場，熟悉市場又懂技術的當老板去了”，世界就是這樣，沒有人能搞定所有問題，所以合作與互補很重要，這就是我為什么一直找志同道合的人的原因。
更多交流請關注我的微信公眾號blacker，言歸正傳，下面這篇文章是劉老師(icoolmedia)寫的，kelly進行編輯和整理
一、譜減法語音降噪基本原理
譜減算法為最早的語音降噪算法之一，它的提出，基於一個簡單的原理：
假設語音中的噪聲只有加性噪聲，只要將帶噪語音譜減去噪聲譜，就可以得到純凈語音幅度。這么做的前提是噪聲信號是平穩的或者緩慢變化的。
得到純凈信號的幅度譜后，可以結合帶噪語音相位（近似帶替純凈語音相位），從而得到近似的純凈語音，語音信號的相位對語音可懂度不敏感。
按上述所示，如果我們設y(n)為受噪聲污染的信號，則y(n)由純凈語音信號x(n)和加性噪聲d(n)組成，即：y(n)=X(n)+d(n)。
其傅里葉變換后表示為：Y(ω)=X(ω)+D(ω)，或寫為：
X(ω) = Y(ω) – D(ω)，如果用功率譜表示可以寫為：
這里
被稱為交叉項，我們假定d(n)具有0均值，並且與x(n)不相關，則交叉項為0，
上述公式簡化為：
或寫為：
二、音樂噪聲和過減因子、譜下限的關系
如果帶噪語音的幅度譜（功率譜也同此理）與估計出來的噪聲譜相減出現負值時，說明對噪聲出現了過估計問題，對這種現象最簡單的處理就是將負值設為0，以保證非負的幅度譜。但是對負值的這種處理，會導致信號幀頻譜的隨機位置上出現小的，獨立的峰值。
轉換到頻域后，這些峰值聽起來就像幀與幀之間頻率隨機變化的多頻音，這種情況在清音段尤其明顯，這種由於半波整流引起的“噪聲”被稱為“音樂噪聲”。從根本上，通常導致音樂噪聲的原因主要有：
1，對譜減算法中的負數部分進行了非線性處理
2，對噪聲譜的估計不准
3，抑制函數（增益函數）具有較大的可變性
減小音樂噪聲的方法是對噪聲譜使用過減技術，同時對譜減后的負值設置一個下限，而不是將它們設為0，其技術形式如下：
其中alpha（大於等於1）為過減因子，它主要影響語音譜的失真程度。Beta（大於0小於1）是譜下限參數，可以控制殘留噪聲的多少以及音樂噪聲的大小。
使用過減因子與譜下限的動機在於：當從帶噪語音譜中減去噪聲譜估計的時候，頻譜中會殘留一些隆起的部分或譜峰，有些譜峰是寬帶的，有些譜峰很窄，看起來像是頻譜上的一個脈沖。通過對噪聲譜的過減處理，我們可以減小寬帶譜峰的幅度，有時還可以將其完全消除。但是僅僅這樣還不夠，因為譜峰周圍可能還存在較深的譜谷。因此需要使用譜下限來“填充”這些譜谷。
在高信噪比中，alpha應取小值；對低信噪比中，alpha建議取大值。Berouti等人做了大量實驗來確定alpha與beta的最優值，在這里我們直接使用就可以了。具體請參考論文：Enhancement of speech corrupted by a acoustic noise。
下面給出譜減算法的Matlab驗證代碼，調用方法為specsub(‘filename.wav’,’outfile.wav’);
1 function specsub(filename,outfile) 2 3 if nargin < 2 4 5 fprintf('Usage: specsub noisyfile.wav outFile.wav \n\n'); 6 7 return; 8 9 end 10 11 12 [x,fs,nbits] = wavread(filename); 13 14 len = floor(20*fs/1000); % Frame size in samples 15 16 if rem(len,2) == 1, len=len+1; end; 17 18 PERC = 50; % window overlap in percent of frame size 19 20 len1 = floor(len*PERC/100); 21 22 len2 = len-len1; 23 24 25 Thres = 3; % VAD threshold in dB SNRseg 26 27 Expnt = 2.0; % power exponent 28 29 beta = 0.002; 30 31 G = 0.9; 32 33 34 win = hamming(len); 35 36 winGain = len2/sum(win); % normalization gain for overlap+add with 50% overlap 37 38 39 % Noise magnitude calculations - assuming that the first 5 frames is noise/silence 40 41 nFFT = 2*2^nextpow2(len); 42 43 noise_mean = zeros(nFFT,1); 44 45 j=1; 46 47 for k = 1:5 48 49 noise_mean = noise_mean+abs(fft(win.*x(j:j+len-1),nFFT)); 50 51 j = j+len; 52 53 end 54 55 noise_mu = noise_mean/5; 56 57 58 %--- allocate memory and initialize various variables 59 60 k = 1; 61 62 img = sqrt(-1); 63 64 x_old = zeros(len1,1); 65 66 Nframes = floor(length(x)/len2)-1; 67 68 xfinal = zeros(Nframes*len2,1); 69 70 71 72 %========================= Start Processing =============================== 73 74 for n = 1:Nframes 75 76 insign = win.*x(k:k+len-1); % Windowing 77 78 spec = fft(insign,nFFT); % compute fourier transform of a frame 79 80 sig = abs(spec); % compute the magnitude 81 82 %save the noisy phase information 83 84 theta = angle(spec); 85 86 SNRseg = 10*log10(norm(sig,2)^2/norm(noise_mu,2)^2); 87 88 if Expnt == 1.0 % 幅度譜 89 90 alpha = berouti1(SNRseg); 91 92 else 93 94 alpha = berouti(SNRseg); % 功率譜 95 96 end 97 98 %&&&&&&&&& 99 100 sub_speech = sig.^Expnt - alpha*noise_mu.^Expnt; 101 102 diffw = sub_speech - beta*noise_mu.^Expnt; % 當純凈信號小於噪聲信號的功率時 103 104 % beta negative components 105 106 z = find(diffw <0); 107 108 if~isempty(z) 109 110 sub_speech(z) = beta*noise_mu(z).^Expnt; % 用估計出來的噪聲信號表示下限值 111 112 end 113 114 % --- implement a simple VAD detector -------------- 115 116 if (SNRseg < Thres) % Update noise spectrum 117 118 noise_temp = G*noise_mu.^Expnt+(1-G)*sig.^Expnt; % 平滑處理噪聲功率譜 119 120 noise_mu = noise_temp.^(1/Expnt); % 新的噪聲幅度譜 121 122 end 123 124 % flipud函數實現矩陣的上下翻轉，是以矩陣的“水平中線”為對稱軸 125 126 %交換上下對稱元素 127 128 sub_speech(nFFT/2+2:nFFT) = flipud(sub_speech(2:nFFT/2)); 129 130 x_phase = (sub_speech.^(1/Expnt)).*(cos(theta)+img*(sin(theta))); 131 132 % take the IFFT 133 134 xi = real(ifft(x_phase)); 135 136 % --- Overlap and add --------------- 137 138 xfinal(k:k+len2-1)=x_old+xi(1:len1); 139 140 x_old = xi(1+len1:len); 141 142 k = k+len2; 143 144 end 145 146 wavwrite(winGain*xfinal,fs,16,outfile); 147 148 function a = berouti1(SNR) 149 150 if SNR >= -5.0 & SNR <= 20 151 152 a = 3-SNR*2/20; 153 154 else 155 156 if SNR < -5.0 157 158 a = 4; 159 160 end 161 162 if SNR > 20 163 164 a = 1; 165 166 end 167 168 end 169 170 function a = berouti(SNR) 171 172 if SNR >= -5.0 & SNR <= 20 173 174 a = 4-SNR*3/20; 175 176 else 177 178 if SNR < -5.0 179 180 a = 5; 181 182 end 183 184 if SNR > 20 185 186 a = 1; 187 188 end 189 190 end 191 192
三、幾種改進的譜減算法
1，非線性譜減
Berouti等人提出的譜減算法，假設了噪聲對所有的頻譜分量都有同等的影響，繼而只用了一個過減因子來減去對噪聲的過估計。現實世界中的噪聲並非如此，這意味着可以用一個頻率相關的減法因子來處理不同類型的噪聲。
2，多帶譜減法
在多帶算法中，將語音頻譜划分為N個互不重疊的子帶，譜減法在每個子帶獨立運行。將語音信號分為多個子帶信號的過程可以通過在時域使用帶通濾波器來進行，或者在頻域使用適當的窗。通常會采用后一種辦法，因為實現起來有更小的運算量。
多帶譜減與非線性譜減的主要區別在於對過減因子的估計。多帶算法針對頻帶估計減法因子，而非線性譜減算法針對每一個頻點，導致頻點上的信噪比可能有很大變化。這種劇烈變化是譜減法中所遇到的語音失真（音樂噪聲）的原因之一。相反，子帶信噪比變化則不會特別劇烈。
3，MMSE譜減算法
上面的方法中，譜減參數alpha和beta通過實驗確定，無論如何都不會是最優的選擇。MMSE譜減法能夠在均方意義下最優地選擇譜減參數。具體請參考論文：A parametic formulation of the generalized spectral subtractor method
4，擴展譜減法
基於自適應維納濾波與譜減原理的結合。維納濾波用於估計噪聲譜，然后從帶噪語音信號中減去該噪聲譜。具體請參考以下兩篇論文：
Extended Spectral Substraction:Description and Preliminary Results.
Extended Spectral Substraction
5，自適應增益平均的譜減
譜減法中導致音樂噪聲的兩個因素在於譜估計的大范圍變化以及增益函數的不同。對於第一個問題，Gustafsson等人建議將分析幀划分為更換小的子幀以得到更低分辨率的頻譜。子幀頻譜通過連續平均以減小頻譜的波動。對於第二個問題Gustafsson等人提出使用自適應指數平均，在時間上對增益函數做平滑。此外，為了避免因使用零相位增益函數導致的非因果濾波問題，Gustafsson等人建議在增益函數中引入線性相位。具體請參考論文：Spectral subtraction using reduced delay convolution and adaptive averaging
6，選擇性譜減法
前面提到的方法對所有語音都做同樣處理。並不區分是濁音段還是清音段。區分濁音與清音的譜減法有：
6.1， 雙頻帶譜減法。通過將帶噪語音能量與某一閾值進行比較，把語音幀分為濁音和清音。對於濁音幀，用算法確定一個截止頻率，在該截止頻率之上，語音被認為是隨機信號。濁音段則通過濾波分為兩個頻帶，一個頻帶位於截止頻率之下（低通濾波后的語音），另外一個頻帶高於截止頻率（高通濾波后的語音）。然后對低通和高通后的語音信號使用不同的算法進行處理。對低通語音部分在短時傅立葉變換的基礎上使用過減算法，對於高通部分以及清音段，使用Thomson的多窗譜估計器取代FFT估計器。主要目的在於減小高頻部分的頻譜值的波動。具體請參考論文：Adaptive two-band spectral subtraction with multi-window spectral estimation
6.2，雙激勵語音模型法，該算法把語音分為兩個獨立的組成部分--濁音分量和清音分量。也就是說，語音由這兩個分量的和來表示（注意不同於將語音分為濁音段和清音段）。濁音分量的分析是基於對基音頻率和諧波幅度的提取。然后從帶噪語音譜中減去濁音譜就得到了清音譜。然后使用一個雙通道系統，基中一個包括改進的維納濾波器，被用於增強清音譜。最終增強的語音由增強后的濁音分量和清音分量求和得到。具體請參考論文：Speech enhancement using the dual excitation speech model
6.3，還有一種基於濁音、清音的譜減算法，在該算法中語音幀首先根據能量和過零率被划分為濁音和清音。然后將帶噪語音譜與銳化函數進行卷積，清音的頻譜就會被銳化（用銳化函數進行鐠銳化的目的在於增加譜對比度，即在抑制譜谷的同時使譜峰更加突出）。具體請參考論文：Spectral subtraction based on phonetic dependency and masking effects
7，基於感知特性的譜減
前面提到的方法，譜減參數要么是通過實驗計算短時信噪比得到，要么是通過最優均方誤差得到，均沒有考慮聽覺系統的特性，該算法的主要目的是使殘余噪聲在聽覺上難以被察覺。利用了人類聽覺系統改進系統的可懂度（即人耳的掩蔽效應）
喜歡系列文章請關注微信公眾號blacker，或掃描下方二維碼：
