# Spectral Subtraction python implementation error

I am trying to implement Spectral subtraction in python. I am new to DSP. I used this article as reference and below code is what I can come up with. Of course I did not get the expected output. I have tried my best to explain the code. Please help me correct the code. Thank you in advance.

def OnSpecSubstract(data, f_size, n_frames):

"""
data is the complete signal
f_size is the size of frames
n_frames is the number of noisy frames
"""

samples = []
phases = []
out = [];

"""
lps is basically the total number of frames that wiil be created, frames will overlap
by 50%
"""
lps = int(len(data)/f_size)+1;

"""
Hanning Window for windowing function
"""
hn_win = np.hanning(f_size)

for i in range(0, (lps*2)):

"""
The 2 line of codes below extract frames of size f_size and append each frame
in list samples. Each frame overlaps by 50% from pervoius frame.
"""
f_loc = int(i*f_size/2);
samples.append( np.asarray( data[ (f_loc):(f_loc+f_size) ] ) )

"""
l_ele is the index of the last element in list samples ie. that frame that was
just appended above.
We perform windowing function on each frame.
Next comes fft of each frame.
Next we save the phase information spererately of each frame in list phases.
In last line we extract the magnitude information of each frame.
"""
l_ele = len(samples)-1;
samples[l_ele] = samples[l_ele] * hn_win[:len( samples[l_ele] )];
samples[l_ele] = fft(samples[l_ele])
phases.append(np.angle(samples[l_ele]))
samples[l_ele] = abs(samples[l_ele])

"""
We check if the current frame has size less than f_size we know that all the
data in signal has been windowed and we break the loop.
"""
if ( len(samples[len(samples)-1]) < f_size ):
break

"""
We assume that first few frames are only noisy frame. Hence, we first initailize
noise by saving the first frame in noise.
"""
noise = np.asarray(samples[0])

"""
Next we add all the noisy frames together. n_frames denotes the number of noisy
frames.
"""
for i in range(1, n_frames):
noise = list(map(lambda x, y : x + y, noise, samples[i]))

"""
Here we take the mean of all the values in noise. We now have the single value for
noise
"""
noise = list(map(lambda x : x/n_frames, noise))
noise = np.mean(noise)
noise *= 2 #The 2 here is bias

"""
We delete the frames that represent noise form samples and phase as we don't
need them.
"""
del samples[:n_frames]
del phases[:n_frames]

"""
Here, we subtract the average noise from each value in each frame.
"""
samples = [ list(map( lambda x : x - noise, i )) for i in samples ]

"""
In the loop below we
First zero off all the negitive amplitudes.
Second, we multiply with each frame their respective phases that was saved earlier.
Third, we inverse fast fourier transform each frame.
"""
for idx, data in enumerate(samples):
samples[idx] = [ 0 if i < 0 else i for i in data ]
samples[idx] = list(map( lambda x, y: x*y, samples[idx], phases[idx] ))
samples[idx] = ifft(samples[idx])

"""
out will be the final reconstructed signal. We first initialize it by saving the
first frame from samples. We then delete the first frame from samples.
"""
out = samples[0]
del samples[0]

"""
Here I try to reconstruct the orignal signal. This is one section where I think there
is definately an issue.
In the below loop what is done is simply the overlapping parts of the frames are
summed up.
"""
for data in samples:

if len(data) == f_size:

"""
dm_pnt is the mid index of data and om_pnt in the index of put that is
f_size/2 places behind the last element.
"""
dm_pnt = int(f_size/2)
om_pnt = int(len(out)-dm_pnt)

out[om_pnt:] = list(map( lambda x, y: x + y, out[om_pnt:], data[:dm_pnt] ))

out = np.concatenate([out, data[dm_pnt:]])

return out

• "I did not get the expected output" so, ok, what did you expect, what did you get instead, how does it differ? What have you considered so far? – Marcus Müller Jul 26 '18 at 7:24

Alright, I was able to figure out a way. I am now getting expected results. However, there is something strange about the output that I don't understand. The result of any Inverst FFT operation is complex numbers ie. real and imaginary part. When we IFFT any signal after processing, to get the actual audio we use only the real part of the IFFT results(correct me if I am wrong). However, here the processed signal is actually the imaginary part. I use the below code to save the imaginary part of the result to a wav file.

1. x = OnSpecSubstract( data, f_size = 400, n_frames = 12 )
2. out = x.imag.astype(np.int16)
3. wavfile.write( 'Output.wav', rate, out)

The real part of the results are very small. The largest real number in the result was around 1.9e-12. Check out this playlist for output.

• The h_samp is noisy data.
• out_b1 is filtered signal with bias = 1.
• out_b2 is filtered signal with bias = 2, and so on.

The amount of musical note increases with bias. I have edited the above code in question to the one that I am using.