对于低频区域语音增强效果很明显,但是无法滤除背景声和其他杂声!下面给出python + matlib版本的实现!

#!/usr/bin/python
from __future__ import division
import numpy as np
import math
from scipy.special import *
import wave
import sys
import matplotlib.pyplot as pltnp.seterr('raise')
infile="audio_8.wav"                               ##输入文件名(读取的音频)
input_file = wave.open(infile, "rb")
params = input_file.getparams()
nchannels, sampwidth, fs, nframes = params[:4]
outfile="noise_reductionnones.wav"                      ##输出文件名(降噪后)
output_file = wave.open(outfile,"wb")
output_file.setnchannels(1)#           设置输出通道
output_file.setsampwidth(2)#           设置信道宽度
output_file.setframerate(2000)#        设置频率S
chunk_size = int(np.floor(60*fs))
#args.initial_noise, args.window_size, args.noise_threshold,
'''
noise_frames:初始化噪音6
Slen:初始化窗口大小0
eta:初始化噪音阈值0.15
'''
def logmmse(x, Srate, noise_frames=6, Slen=0, eta=0.15, saved_params=None):
#################################初始化变量########################################if Slen == 0:Slen = int(math.floor(0.02 * Srate))#40################################################保证Slen为偶数if Slen % 2 == 1:Slen = Slen + 1PERC = 50#                                      窗口重叠占帧大小的百分比#math.floor向下取整5.5->5len1 = math.floor(Slen * PERC / 100)#20         样本中的帧大小len2 = int(Slen - len1)#20                      样本中的更新率#W(n,α ) = (1 -α ) - α cos(2*PI*n/(N-1)),0≦n≦N-1win = np.hanning(Slen)#                         定义窗口win = win * len2 / np.sum(win)#                 normalise窗口用于等电平输出nFFT = 2 * Slen#分配内存并初始化各种变量x_old = np.zeros(len1)Xk_prev = np.zeros(len1)Nframes = int(math.floor(len(x) / len2) - math.floor(Slen / len2))xfinal = np.zeros(Nframes * len2)noise_mean = np.zeros(nFFT)for j in range(0, Slen*noise_frames, Slen):noise_mean = noise_mean + abs(np.fft.fft(win * x[j:j + Slen], nFFT, axis=0))noise_mu2 = noise_mean / noise_frames ** 2'''if saved_params == None:noise_mean = np.zeros(nFFT)for j in range(0, Slen*noise_frames, Slen):noise_mean = noise_mean + abs(np.fft.fft(win * x[j:j + Slen], nFFT, axis=0))noise_mu2 = noise_mean / noise_frames ** 2else:noise_mu2 = saved_params['noise_mu2']Xk_prev = saved_params['Xk_prev']x_old = saved_params['x_old']'''print('saved_params',saved_params)
###############################开始处理##############################################################aa = 0.98mu = 0.98eta = 0.15ksi_min = 10 ** (-25 / 10)cnt=0for k in range(0, Nframes*len2, len2):insign = win * x[k:k + Slen]spec = np.fft.fft(insign, nFFT, axis=0)cnt +=1if k==0:print('spec',spec[0])sig = np.absolute(spec)#                           计算幅度sig2 = sig ** 2gammak = np.minimum(sig2 / noise_mu2, 40)if Xk_prev.all() == 0:ksi = aa + (1 - aa) * np.maximum(gammak - 1, 0)# 则在SNR后限制,以避免溢出else:ksi = aa * Xk_prev / noise_mu2 + (1 - aa) * np.maximum(gammak - 1, 0)  #a先验SNRksi = np.maximum(ksi_min, ksi)log_sigma_k = gammak * ksi/(1 + ksi) - np.log(1 + ksi)vad_decision = np.sum(log_sigma_k)/Slenif (vad_decision < eta):#噪声帧仅发现noise_mu2 = mu * noise_mu2 + (1 - mu) * sig2##############################结束##############################################################A = ksi / (1 + ksi)                                                        #Log-MMSE估计器vk = A * gammakei_vk = 0.5 * expn(1, vk)hw = A * np.exp(ei_vk)sig = sig * hwXk_prev = sig ** 2xi_w = np.fft.ifft(hw * spec, nFFT, axis=0)xi_w = np.real(xi_w)xfinal[k:k + len2] = x_old + xi_w[0:len1]x_old = xi_w[len1:Slen]print('cnt',cnt)return xfinal, {'noise_mu2': noise_mu2, 'Xk_prev': Xk_prev, 'x_old': x_old}
frames_read = 0
while (frames_read < nframes):frames = nframes - frames_read if frames_read + chunk_size > nframes else chunk_size# signal = input_file.read_frames(frames)signal = input_file.readframes(frames)signal = np.fromstring(signal, dtype='int16')signal=signal/15000frames_read = frames_read + framesoutput, saved_params = logmmse(signal, fs)output = np.array(output*np.iinfo(np.int16).max, dtype=np.int16)output_file.writeframes(output)
input_file.close()
output_file.close()
function logmmse(filename,outfile)%
%  Implements the logMMSE algorithm [1].
%
%  Usage:  logmmse(noisyFile, outputFile)
%
%         infile - noisy speech file in .wav format
%         outputFile - enhanced output file in .wav format
%
%
%  Example call:  logmmse('sp04_babble_sn10.wav','out_log.wav');
%
%  References:
%   [1] Ephraim, Y. and Malah, D. (1985). Speech enhancement using a minimum
%       mean-square error log-spectral amplitude estimator. IEEE Trans. Acoust.,
%       Speech, Signal Process., ASSP-23(2), 443-445.
%
% Authors: Philipos C. Loizou
%
% Copyright (c) 2006 by Philipos C. Loizou
% $Revision: 0.0 $  $Date: 10/09/2006 $
%-------------------------------------------------------------------------if nargin<2fprintf('Usage: logmmse(noisyfile.wav,outFile.wav) \n\n');return;
end[x, Srate, bits]= wavread( filename);   %nsdata is a column vector% =============== Initialize variables ===============len=floor(20*Srate/1000); % Frame size in samples
if rem(len,2)==1, len=len+1; end;
PERC=50; % window overlap in percent of frame size
len1=floor(len*PERC/100);
len2=len-len1;win=hamming(len);  % define window% Noise magnitude calculations - assuming that the first 6 frames is
% noise/silence nFFT=2*len;
noise_mean=zeros(nFFT,1);
j=1;
for m=1:2noise_mean=noise_mean+abs(fft(win.*x(j:j+len-1),nFFT));j=j+len;
end
noise_mu=noise_mean/6;
noise_mu2=noise_mu.^2;%--- allocate memory and initialize various variablesx_old=zeros(len1,1);
Nframes=floor(length(x)/len2)-floor(len/len2);
xfinal=zeros(Nframes*len2,1);%===============================  Start Processing =======================================================
%
k=1;
aa=0.98;
mu=0.98;
eta=0.15; ksi_min=10^(-25/10);for n=1:Nframesinsign=win.*x(k:k+len-1);spec=fft(insign,nFFT);sig=abs(spec); % compute the magnitudesig2=sig.^2;gammak=min(sig2./noise_mu2,40);  % limit post SNR to avoid overflowsif n==1ksi=aa+(1-aa)*max(gammak-1,0);elseksi=aa*Xk_prev./noise_mu2 + (1-aa)*max(gammak-1,0);     % a priori SNRksi=max(ksi_min,ksi);  % limit ksi to -25 dBendlog_sigma_k= gammak.* ksi./ (1+ ksi)- log(1+ ksi);    vad_decision= sum(log_sigma_k)/ len;    if (vad_decision< eta) % noise only frame foundnoise_mu2= mu* noise_mu2+ (1- mu)* sig2;end% ===end of vad===A=ksi./(1+ksi);  % Log-MMSE estimatorvk=A.*gammak;ei_vk=0.5*expint(vk);hw=A.*exp(ei_vk);sig=sig.*hw;Xk_prev=sig.^2;xi_w= ifft( hw .* spec,nFFT);xi_w= real( xi_w);xfinal(k:k+ len2-1)= x_old+ xi_w(1:len1);x_old= xi_w(len1+ 1: len);k=k+len2;endwavwrite(xfinal,Srate,16,outfile);

logmmse降噪算法相关推荐

  1. 可复现的图像降噪算法总结——超赞整理

    点击上方"小白学视觉",选择加"星标"或"置顶" 重磅干货,第一时间送达 本文转自:AI算法与图像处理 图像降噪,是最简单也是最基础的图像处 ...

  2. 图像降噪算法——低秩聚类:WNNM算法

    图像降噪算法--低秩聚类:WNNM算法 图像降噪算法--低秩聚类:WNNM算法 1. 基本原理 2. matlab代码 3. 结论 图像降噪算法--低秩聚类:WNNM算法 同样是为了完善自己知识版图的 ...

  3. 图像降噪算法——稀疏表达:K-SVD算法

    图像降噪算法--稀疏表达:K-SVD算法 图像降噪算法--稀疏表达:K-SVD算法 1. 基本原理 2. python代码 3. 结论 图像降噪算法--稀疏表达:K-SVD算法 为了完善下自己降噪算法 ...

  4. 图像降噪算法——DnCNN / FFDNet / CBDNet / RIDNet / PMRID / SID

    图像降噪算法--DnCNN / FFDNet / CBDNet / RIDNet / PMRID / SID 图像降噪算法--DnCNN / FFDNet / CBDNet / RIDNet / PM ...

  5. 图像降噪算法——Variance Stabilizing Transform / Generalization Anscombe Transform算法

    图像降噪算法--Variance Stabilizing Transform / Generalization Anscombe Transform算法 图像降噪算法--Variance Stabil ...

  6. 图像降噪算法——时域降噪算法

    图像降噪算法--时域降噪算法 图像降噪算法--时域降噪算法 1.<MeshFLow Video Denoising> 2. <Real-Time Video Denoising On ...

  7. 图像降噪算法——小波硬阈值滤波(下)

    图像降噪算法--小波硬阈值滤波(下) 图像降噪算法--小波硬阈值滤波(下) 1. 基本原理 2. C++代码实现 3. 结论 图像降噪算法--小波硬阈值滤波(下) 1. 基本原理 关于离散小波变换的原 ...

  8. 图像降噪算法——小波硬阈值滤波(上)

    图像降噪算法--小波硬阈值滤波(上) 图像降噪算法--小波硬阈值滤波(上) 1. 多分辨率展开 2. 尺度函数 3. 小波函数 4. 小波级数展开 5. 离散小波变换 6. 快速小波变换 7. 图像小 ...

  9. 图像降噪算法——维纳滤波

    图像降噪算法--维纳滤波 图像降噪算法--维纳滤波 1. 基本原理 2. C++代码实现 3. 结论 图像降噪算法--维纳滤波 维纳滤波是在频域中处理图像的一种算法,是一种非常经典的图像增强算法,不仅 ...

  10. 图像降噪算法——高斯低通滤波

    图像降噪算法--高斯低通滤波 图像降噪算法--高斯低通滤波 1. 基本原理 2. C++代码实现 3. 结论 图像降噪算法--高斯低通滤波 1. 基本原理 通过离散傅里叶变换对图像进行滤波流程作非常简 ...

最新文章

  1. 梯度下降优化算法综述与PyTorch实现源码剖析
  2. sdut 2805(最小生成树)
  3. DB2 导出数据文件
  4. vs2015编译 pybind 动态库
  5. 解决IntelliJ IDEA 2019.3.5 启动无反应
  6. 网络编程练习 -- NSURLConnection -- get/post请求
  7. C和C++混合编程(__cplusplus使用)
  8. textmate bundle for jquery
  9. java 类数组_Java常用类-字符串、日期类、算法及数组工具类等
  10. 图像 super-resolution restruction 的各种主流实现方式
  11. 旧版sai笔刷_最详细的SAI笔刷设置教程,非常全面详细!
  12. java静态链表_Java数据结构——静态链表实现
  13. 服务器开机自检提示信息,hp服务器开机自检报错提示
  14. Oracle的SQL注入
  15. bigdecimal负数变正数_Java中BigDecimal的8种舍入模式
  16. javaScript深度克隆
  17. coredata 详解
  18. 【异常记录】C# 连接数据库错误异常解决 err:40
  19. 百度服务器自动重启,百度云服务器重启的两种方法介绍
  20. 广东高中生多少人_2020年广东高考报名人数统计有多少人

热门文章

  1. set+线段树 Codeforces Round #305 (Div. 2) D. Mike and Feet
  2. hdu 1284 钱币兑换问题 (递推 || DP || 母函数)
  3. 驱动编译的时候注意编译工程选项
  4. 用户体验设计师、UI 设计师和交互设计师之间的区别,如何挑选图书?
  5. [2018.08.07 T1] 签到?
  6. 计算机硬盘分区知识简介、Linux企业级分区方案建议
  7. delete postman 传参_Postman高级应用——串行传参和动态传参详解
  8. 合肥青少年信息学计算机竞赛试题,合肥市第二十九届青少年信息学奥林匹克竞赛(小学组)试题及部分答案...
  9. matlab在二值图像上画曲线_数字图像处理:Image Printing Program Based on Halftoning
  10. html头部打开页面为兼容模式,Web页面因为兼容模式产生的奇怪问题解答