How to extract Audio Wave from a mixture of Signal using Scipy – Python?
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:
Attention geek! Strengthen your foundations with the Python Programming Foundation Course and learn the basics.
To begin with, your interview preparations Enhance your Data Structures concepts with the Python DS Course. And to begin with your Machine Learning Journey, join the Machine Learning - Basic Level Course
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 signalfrom scipy import signalfrom scipy.fft import fftimport numpy as npimport matplotlib.pyplot as plt # Number of samplesN_sample = 512 # Sampling frequencyfs = 10000 # inter sample time = 0.001s = 1kHz samplingdt = 1/fs # time vectort = np.arange(0, N_sample)*dt # Create signal vector that is the sum of 50 Hz, 1023 Hz, and 1735 HzSignal = 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 signalSignal = Signal+np.random.normal(0, .1, Signal.shape) # Part A: Estimation of Length and Window# Select design Specification# Lower stopband frequency in Hzfstop_L = 500 # Lower passband frequency in HZfpass_L = 800 # Upper stopband frequency in Hzfstop_U = 1500 # Upper passband frequency in HZfpass_U = 1200 # Calculations# Normalized lower transition band w.r.t. fsdel_f1 = abs(fpass_L-fstop_L)/fs # Normalized upper transition band w.r.t. fsdel_f2 = abs(fpass_U-fstop_U)/fs # Filter length using selected window based# on Normalized lower transition bandN1 = 3.3/del_f1 # Filter length using selected window based# on Normalized upper transition bandN2 = 3.3/del_f2print('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 integerN = int(np.ceil(max(N1, N2)))if(N % 2 == 0): N = N+1print('Selected filter length :', N) # Calculate lower and uper cut-off frequencies# Lower cut-off frequency in HzfL = (fstop_L+fpass_L)/2 # Upper cut-off frequency in HzfU = (fstop_U+fpass_U)/2 # Normalized Lower cut-off frequency in (w/pi) radwL = 2*fL/fs # Normalized upper cut-off frequency in (w/pi) radwU = 2*fU/fs # Cutoff frequency arraycutoff = [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 windowb_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 windoww, h_ham = signal.freqz(b_ham, a=1) # Calculate Magnitude in dB# Calculate magnitude in decibelsh_dB_ham = 20*np.log10(abs(h_ham)) a = [1] # Filter the noisy signal by designed filter# using signal.filtfiltfiltOut = signal.filtfilt(b_ham, a, Signal) # Plot filter magnitude and phase responses using# subplot. Digital frequency w converted in analog# frequencyfig = plt.figure(figsize=(12, 18)) # Original signalsub1 = 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 Plotsub2 = 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 windowsub3.set_ylabel('Magnitude ')sub3.set_xlabel('Time')sub3.set_title('Filtered output of Band pass Filter', fontsize=20) # Show spectrum of noisy input signalSigf = fft(Signal) # Compute FFT of noisy signalsub4 = 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 signalOutf = fft(filtOut) # Compute FFT of filtered signalsub5 = 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() |
Output:
