Example of Python Spectral Subtraction Speech Noise Reduction

#!/usr/bin/env python
import numpy as np
import wave
import nextpow2
import math
 
 
# Open a WAV file
f = ("")
# Read format information
# (nchannels, sampwidth, framerate, nframes, comptype, compname)
params = ()
nchannels, sampwidth, framerate, nframes = params[:4]
fs = framerate
# Read waveform data
str_data = (nframes)
()
# Convert waveform data to an array
x = (str_data, dtype=)
# Calculated parameters
len_ = 20 * fs // 1000
PERC = 50
len1 = len_ * PERC // 100
len2 = len_ - len1
# Set default parameters
Thres = 3
Expnt = 2.0
beta = 0.002
G = 0.9
# Initialize the Hamming window
win = (len_)
# normalization gain for overlap+add with 50% overlap
winGain = len2 / sum(win)
 
 
# Noise magnitude calculations - assuming that the first 5 frames is noise/silence
nFFT = 2 * 2 ** (nextpow2.nextpow2(len_))
noise_mean = (nFFT)
 
 
j = 0
for k in range(1, 6):
  noise_mean = noise_mean + abs((win * x[j:j + len_], nFFT))
  j = j + len_
noise_mu = noise_mean / 5
 
 
# --- allocate memory and initialize various variables
k = 1
img = 1j
x_old = (len1)
Nframes = len(x) // len2 - 1
xfinal = (Nframes * len2)
 
 
# =========================  Start Processing  ===============================
for n in range(0, Nframes):
  # Windowing
  insign = win * x[k-1:k + len_ - 1]
  # compute fourier transform of a frame
  spec = (insign, nFFT)
  # compute the magnitude
  sig = abs(spec)
 
 
  # save the noisy phase information
  theta = (spec)
  SNRseg = 10 * np.log10((sig, 2) ** 2 / (noise_mu, 2) ** 2)
 
 
 
 
  def berouti(SNR):
    if -5.0 <= SNR <= 20.0:
      a = 4 - SNR * 3 / 20
    else:
      if SNR < -5.0:
        a = 5
      if SNR > 20:
        a = 1
    return a
 
 
 
 
  def berouti1(SNR):
    if -5.0 <= SNR <= 20.0:
      a = 3 - SNR * 2 / 20
    else:
      if SNR < -5.0:
        a = 4
      if SNR > 20:
        a = 1
    return a
 
 
  if Expnt == 1.0: # Amplitude spectrum
    alpha = berouti1(SNRseg)
  else: # Power spectrum
    alpha = berouti(SNRseg)
  #############
  sub_speech = sig ** Expnt - alpha * noise_mu ** Expnt;
  # When the pure signal is less than the power of the noise signal
  diffw = sub_speech - beta * noise_mu ** Expnt
  # beta negative components
 
 
  def find_index(x_list):
    index_list = []
    for i in range(len(x_list)):
      if x_list[i] < 0:
        index_list.append(i)
    return index_list
 
 
  z = find_index(diffw)
  if len(z) > 0:
    # Express the lower limit in terms of the estimated noise signal
    sub_speech[z] = beta * noise_mu[z] ** Expnt
    # --- implement a simple VAD detector --------------
    if SNRseg < Thres: # Update noise spectrum
      noise_temp = G * noise_mu ** Expnt + (1 - G) * sig ** Expnt # Smoothed processing noise power spectrum
      noise_mu = noise_temp ** (1 / Expnt) # New noise amplitude spectrum
    # The flipud function flips the matrix up and down, using the matrix's "horizontal centerline" as the axis of symmetry.
    # Swap up and down symmetrical elements
    sub_speech[nFFT // 2 + 1:nFFT] = (sub_speech[1:nFFT // 2])
    x_phase = (sub_speech ** (1 / Expnt)) * (([(x) for x in theta]) + img * (([(x) for x in theta])))
    # take the IFFT
 
 
    xi = (x_phase).real
    # --- Overlap and add ---------------
    xfinal[k-1:k + len2 - 1] = x_old + xi[0:len1]
    x_old = xi[0 + len1:len_]
    k = k + len2
# Save the document
wf = ('', 'wb')
# Setting parameters
(params)
# set the wave file .tostring() converts array to data
wave_data = (winGain * xfinal).astype()
(wave_data.tostring())
()