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 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 upper 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 windowsub2.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:
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 upper 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 windowsub2.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() |
