How to extract Audio Wave from a mixture of Signal using Scipy - Python?

How to extract Audio Wave from a mixture of Signal using Scipy – Python?

Last Updated : 24 Feb, 2021

Prerequisites: Scipy

Spectral Analysis refers to analyzing of the frequency spectrum/response of the waves. This article as the title suggests deals with extracting audio wave from a mixture of signals and what exactly goes into the process can be explained as:

Consider we have 3 mixed Audio Signals having frequency of 50Hz,1023Hz & 1735Hz respectively. Apart from these signals we will be also implementing noise to the signal beforehand. The spectral analysis will be done via using a filter so that we can separate out the signals. On requirement, we can tweak our signals according to the frequency of the signal we want to extract.

Approach

Import modules
Specify conditions such as number of samples, sampling frequency, inner sample time & creating our mixed audio wave
Add noise to the audio signal
Estimate of Filter Window & Computing Cutoff Frequency
Create a filter to filter out noise
Plot the Noisy Signal, Frequency Response of Filter, Extracted Audio Wave, Frequency Spectrum of Mixed Audio Signal, Frequency Spectrum of our extracted Audio Signal
Display plot

Program:

Python3

# Original, high sample rate signal
# Let us imagine this is like our analog signal
from scipy import signal
from scipy.fft import fft
import numpy as np
import matplotlib.pyplot as plt
  
# Number of samples
N_sample = 512
  
# Sampling frequency
fs = 10000
  
# inter sample time = 0.001s = 1kHz sampling
dt = 1/fs
  
# time vector
t = np.arange(0, N_sample)*dt
  
# Create signal vector that is the sum of 50 Hz, 1023 Hz, and 1735 Hz
Signal = np.sin(2*np.pi*50*t) + np.sin(2*np.pi*1023*t)+np.sin(2*np.pi*1735*t)
  
# Add random noise to the signal
Signal = Signal+np.random.normal(0, .1, Signal.shape)
  
# Part A: Estimation of Length and Window
# Select design  Specification
# Lower stopband frequency in Hz
fstop_L = 500
  
# Lower passband frequency in HZ
fpass_L = 800
  
# Upper stopband frequency in Hz
fstop_U = 1500
  
# Upper passband frequency in HZ
fpass_U = 1200
  
# Calculations
# Normalized lower transition band w.r.t. fs
del_f1 = abs(fpass_L-fstop_L)/fs
  
# Normalized upper transition band w.r.t. fs
del_f2 = abs(fpass_U-fstop_U)/fs
  
# Filter length using selected window based
# on Normalized lower transition band
N1 = 3.3/del_f1
  
# Filter length using selected window based
# on Normalized upper transition band
N2 = 3.3/del_f2
print('Filter length based on lower transition band:', N1)
print('Filter length based on upper transition band:', N2)
  
# Select length as the maximum of the N1 and N2
# and if it is even, make it next higher integer
N = int(np.ceil(max(N1, N2)))
if(N % 2 == 0):
    N = N+1
print('Selected filter length :', N)
  
# Calculate lower and uper cut-off frequencies
# Lower cut-off frequency in Hz
fL = (fstop_L+fpass_L)/2
  
# Upper cut-off frequency in Hz
fU = (fstop_U+fpass_U)/2
  
# Normalized Lower cut-off frequency in (w/pi) rad
wL = 2*fL/fs
  
# Normalized upper cut-off frequency in (w/pi) rad
wU = 2*fU/fs
  
# Cutoff frequency array
cutoff = [wL, wU]
  
# Since the given specification of Stopband attenuation = 50 dB
# and Passband ripple = 0.05 dB, atleast satisfy with
# Hamming window, we have to choose it.
  
# Determine Filter coefficients
# Call filter design function using Hamming window
b_ham = signal.firwin(N, cutoff, window="hamming", pass_zero="bandpass")
  
# Determine Frequency response of the filters
# Calculate response h at specified frequency
# points w for Hamming window
w, h_ham = signal.freqz(b_ham, a=1)
  
# Calculate Magnitude in dB
# Calculate magnitude in decibels
h_dB_ham = 20*np.log10(abs(h_ham))
  
a = [1]
  
# Filter the noisy signal by designed filter
# using signal.filtfilt
filtOut = signal.filtfilt(b_ham, a, Signal)
  
  
# Plot filter magnitude and phase responses using
# subplot. Digital frequency w converted in analog
# frequency
fig = plt.figure(figsize=(12, 18))
  
# Original signal
sub1 = plt.subplot(5, 1, 1)
sub1.plot(t[0:200], Signal[0:200])
sub1.set_ylabel('Amplitude')
sub1.set_xlabel('Time')
sub1.set_title('Noisy signal', fontsize=20)
  
# Magnitude response Plot
sub2 = plt.subplot(5, 1, 2)
sub2.plot(w*fs/(2*np.pi), h_dB_ham, 'r', label='Bandpass filter',
          linewidth='2')  # Plot for magnitude response window
  
sub2.set_ylabel('Magnitude (db)')
sub2.set_xlabel('Frequency in Hz')
sub2.set_title('Frequency response of Bandpass Filter', fontsize=20)
sub2.axis = ([0,  fs/2,  -110,  5])
sub2.grid()
  
sub3 = plt.subplot(5, 1, 3)
sub3.plot(t[0:200], filtOut[0:200], 'g', label='Filtered signal',
          linewidth='2')  # Plot for magnitude response window
sub3.set_ylabel('Magnitude ')
sub3.set_xlabel('Time')
sub3.set_title('Filtered output of Band pass Filter', fontsize=20)
  
  
# Show spectrum of noisy input signal
Sigf = fft(Signal)  # Compute FFT of noisy signal
sub4 = plt.subplot(5, 1, 4)
xf = np.linspace(0.0, 1.0/(2.0*dt), (N_sample-1)//2)
sub4.plot(xf, 2.0/N_sample * np.abs(Sigf[0:(N_sample-1)//2]))
sub4.set_ylabel('Magnitude')
sub4.set_xlabel('Frequency in Hz')
sub4.set_title('Frequency Spectrum of Original Signal', fontsize=20)
sub4.grid()
  
# Show spectrum of filtered output signal
Outf = fft(filtOut)  # Compute FFT of filtered signal
sub5 = plt.subplot(5, 1, 5)
xf = np.linspace(0.0, 1.0/(2.0*dt), (N_sample-1)//2)
sub5.plot(xf, 2.0/N_sample * np.abs(Outf[0:(N_sample-1)//2]))
sub5.set_ylabel('Magnitude')
sub5.set_xlabel('Frequency in Hz')
sub5.set_title('Frequency Spectrum of Filtered Signal', fontsize=20)
sub5.grid()
  
fig.tight_layout()
plt.show()