2

I have 2D array of numpy.int16 (44100Hz) audio samples (1376 of them). Each sample has this format (example values):

[ -4 4 -5 -10 -5 -6 -11 -4 -9 -7 -10 1 -4 -8 -9 -8 -4 -13 -14 -11 -12 -4 -14 -13 -9 -2 -2 -16 -5 -5 -4 3 -5 -4 -8 -11 -10 -12 -16 -7 -8 -14 -14 -14 -16 -17 -8 -13 -9 -6 -9 -6 -9 -8 -12 -1 -4 -8 -2 -2 -2 -8 -8 2 1 -8 3 2 0 -6 0 9 0 -2 0 -1 0 -3 -1 1 -2 -2 0 -6 -1 -2 -3 5 -3 1 -1 -5 -3 0 -3 0 -3 -3 -2 -1 -2 1 -4 1 -6 -3 -2 -4 1 1 -1 -6 2 -1 -2 -5 -6 -5 -6 -2 -4 -1 0 -3 -6 -4 -5 -3]

The entire array is called sample_list (yes, very surprising)

I am attempting to perform a bandpass filter on the sample set in the frequency domain, and then convert it back into the time domain format above. This code below converts the samples perfectly from time domain to frequency domain and then back to time domain, although I still need to apply my bandpass filter to it:

import numpy as np

for samp in sample_list:

    time = np.linspace(0,1,44100) # currently not using

    float_samp = np.float32(samp)
    fft_spectrum = np.fft.rfft(float_samp)

    freqs = np.fft.rfftfreq(len(fft_spectrum), d=time[1]-time[0])  # currently not using

    ### need bandpass filter here (brickwall is fine) - 5000Hz to 8000Hz

    time_domain = np.fft.irfft(fft_spectrum)
    converted = np.int16(time_domain)

Once filtered, I plan to run my filtered samples through the code found on this page: Python frequency detection to find the fundamental frequency of the tone. The frequency detection code currently works just fine on my unfiltered samples - although it will give me the frequency of the loudest tone, not necessarily the frequency of the tone i'm looking for... thus my need to filter.

A brick-wall filter is fine because I'm looking for the timing and presence of my test tones, and not necessarily the quality.I'm also concerned about extra processing cycles if they're not needed.

I cannot use Scipy and its fft/filtering library, unfortunately, because I am running the code on Android and Scipy is not available for the code platform I'm using (Kivy). I have to use strictly the Numpy library.

Any help/direction is greatly appreciated.

Community
  • 1
  • 1
A. Wylie - DigiCoyote
  • 67
  • 1
  • 13

1 Answers1

0

If a brick wall filter is acceptable, you can just use,

bw_filter = np.zeros(freqs.shape, dtype='float32')
f_0 = 0.5*(8000 + 5000)
df_0 = 0.5*(8000-5000)
bw_filter[np.abs(freqs - f_0) < df_0] = 1.0

fft_spectrum *= bw_filter

I would have been better to use, for instance scipy.signal.butter, but since you don't have access to scipy, reimplementing it might take time.

rth
  • 10,680
  • 7
  • 53
  • 77
  • 2
    This is not a terribly good brick-wall filter because of the Gibbs Phenomenon. - [http://en.wikipedia.org/wiki/Gibbs_phenomenon]. You're much better off doing your filtering in the time domain. – marko May 01 '15 at 13:11
  • @rth: the shape of freqs is 33, and the shape of fft_spectrum is 65, which threw a _"ValueError: operands could not be broadcast together with shapes (65,) (33,) (65,)"_. I ran this as a test to see it wide open: `bw_filter = np.zeros(65, dtype='float32') f_0 = 0.5*(20000 + 0) df_0 = 0.5*(20000 - 0) bw_filter[np.abs(freqs - f_0) < df_0] = 1.0 fft_spectrum *= bw_filter` It looks like it did process my audio, but if I've put 20000Hz and 0Hz as my band pass, technically it shouldn't change the signal, yes? – A. Wylie - DigiCoyote May 02 '15 at 01:10
  • @marko: Given the choice, I would much rather just process it in the time domain, as will take far less processing power that way instead of converting it into the frequency domain and back again, I'm assuming. What would a bandpass in the time domain look like, code wise, without using scipy? – A. Wylie - DigiCoyote May 02 '15 at 01:17
  • I agree with @marco comment, this is not a great way of doing it. My knowledge of how to do this better pretty much stops at `from scipy.signal import something` and since you don't have access to `scipy` you would have re-implement the corresponding functions or find another library that includes them. – rth May 02 '15 at 07:52
  • @digicoyote: You'd probably do this with an IIR filter. Applying the filter to each sample is convolution - and boils down to a weight sum of samples. See http://en.wikipedia.org/wiki/Infinite_impulse_response. Calculating the filter coefficients is the only difficult part. Matlab is how most people do it, but it's entirely possible to do from first principles. – marko May 02 '15 at 09:45
  • @marko: I can feel that you're itching to post a working python script here. Don't let me stop you from living the dream. ;) – A. Wylie - DigiCoyote May 02 '15 at 19:53
  • @marko: I've gone another direction, looking for magnitude changes in specific bins instead of applying a band-pass filter. I think I should still do a window on the time-domain. – A. Wylie - DigiCoyote May 05 '15 at 07:17