2
$\begingroup$

I'm writing a script for sound-activated recording in Python using pyaudio. I want to trigger a 5s recording after a sound that is above a prespecified volume and frequency. I've managed to get the volume part working but don't know how to specify the minimum trigger frequency (I'd like it to trigger at frequencies above 10kHz, for example):

import pyaudio
import wave
from array import array
import time
 
FORMAT=pyaudio.paInt16
CHANNELS=1
RATE=44100
CHUNK=1024
RECORD_SECONDS=5

audio=pyaudio.PyAudio() 

stream=audio.open(format=FORMAT,channels=CHANNELS, 
                  rate=RATE,
                  input=True,
                  frames_per_buffer=CHUNK)

nighttime=True

while nighttime:
     data=stream.read(CHUNK)
     data_chunk=array('h',data)
     vol=max(data_chunk)
     if(vol>=3000):
         print("recording triggered")
         frames=[]
         for i in range(0, int(RATE / CHUNK * RECORD_SECONDS)):
             data = stream.read(CHUNK)
             frames.append(data)
         print("recording saved")
         # write to file
         words = ["RECORDING-", time.strftime("%Y%m%d-%H%M%S"), ".wav"]
         FILE_NAME= "".join(words) 
         wavfile=wave.open(FILE_NAME,'wb')
         wavfile.setnchannels(CHANNELS)
         wavfile.setsampwidth(audio.get_sample_size(FORMAT))
         wavfile.setframerate(RATE)
         wavfile.writeframes(b''.join(frames))
         wavfile.close()
     # check if still nighttime
     nighttime=True 
 
 stream.stop_stream()
 stream.close()
 audio.terminate()

I'd like to add to the line if(vol>=3000): something like if(vol>=3000 and frequency>10000): but I don't know how to set up frequency. How to do this?

$\endgroup$
1
$\begingroup$

You need to perform FFT on data_chunk, then take the magnitude of the resulting complex array. The result will be a vector with the magnitudes of all frequency components in data_chunk. This is how you would go about doing that:

from scipy.fft import fft
import numpy as np

Y = fft(data_chunk)
Y = Y[0:round(len(Y)/2)]
YMag = np.abs(Y)

Now YMag has the magnitudes of frequencies ranging from 0 Hz, up to $F_s/2$ Hz where $F_s$ is the sampling frequency, i.e. how many samples you receive each second. Now to find a specific frequency in YMag, let's say $f$ you would use the following formula $n = 2*f*N/F_s$, where N is the number of elements in YMag and n is the position of the $f$ frequency component in YMag.

$\endgroup$
5
  • $\begingroup$ Amazing, thanks! I'll try this out later. I think I follow the formula at the end. How would I write the code for if(frequency>10000)? $\endgroup$
    – Thomas
    Jul 14 at 12:51
  • $\begingroup$ I would check whether the magnitude of any of the components higher than the index $n$ have a value above a certain threshold. Something like: if(max(YMag[n:]) > THRESHOLD) $\endgroup$
    – Shivalnu
    Jul 14 at 13:00
  • $\begingroup$ I'm not quite clear on how to pick an index number for n. Would this not work: if(max(YMag) > THRESHOLD)? $\endgroup$
    – Thomas
    Jul 14 at 21:04
  • $\begingroup$ No, this would not work as YMag is an array. You do not pick a number for $n$, you pick a value for $f$ (10000 Hz in your case) and you calculate $n$ from the formula I've written in my answer. $\endgroup$
    – Shivalnu
    Jul 15 at 7:38
  • $\begingroup$ OK I'm struggling with this. n = 2 * 10000 * len(YMag)/44100 doesn't seem to make sense. Could you write the code to get this to work in the if statement that triggers the recording? (if(vol>=3000 and frequency>10000):) $\endgroup$
    – Thomas
    Jul 15 at 8:51
1
$\begingroup$

What you need is to filter out lower frequencies: You do that with a high-pass filter.

Since you don't really need that filter from "very much blocking" to "completely letting everything through" in a very narrow band, that filter will be pretty short (meaning it only needs very few multiplications per sample to compute its output.

I've install pyfda to do a quick and dirty design¹, using the same sample rate as you, a cutoff frequency of 12 kHz, and a transition width of 4 kHz:

Magnitude frequency response

In terms of "how does this filter actually look like":

b = [0.04889, -0.1332, 0.2142, -0.2142, 0.1332, -0.04889]
a = [1,        1.217,  1.996,   1.56,   0.9704,  0.3965 ]

That's it. That's the filter.

You apply it like this: ATTENTION This code was written free-hand and not even run through Python to check for syntactically illegal typos. In other words, you'll have to test it.

from scipy import signal
import numpy
THRESHOLD = 3000
# … the rest of your code up to the while-loop:


"""
Your filter: these are the coefficients for the feedforward part (b) and feedback part (a) of an IIR filter. I got these from pyfda.
"""
b = [0.04889, -0.1332, 0.2142, -0.2142, 0.1332, -0.04889]
a = [1,        1.217,  1.996,   1.56,   0.9704,  0.3965 ]

"""
This is yet another filter, a low-pass filter used for smoothing.
Design specs were fs=44100 Hz, f_bp = 60 Hz, f_sb = 100 Hz,
A_bp = 3 dB, A_sb = 30 db, it's also an elliptic IIR.

The 60 Hz cutoff should get rid of anything that's shorter than 
ca 1/60 of a second – we don't care about things shorter than that,
was my assumption.
"""
b_lpf = [0.0005067106574837601, -0.000506612969002681, -0.0005066129690026811, 0.00050671065748376]
a_lpf = [1.0                  , -2.9949235962214953,    2.9899182657208225, -0.9949944741223649]

"""
Initialize the filter state so that it starts "eingeschwungen" (in steady state)
This is optional, but without it, chances are high you always trigger when you first enter the while loop.
"""
filter_state = signal.lfilter_zi(b, a)
filter_state_lpf = signal.lfilter_zi(b_lpf, a_lpf)

left_to_record = 0
wavfile = None
leftover_data = []

while nighttime:
    data=stream.read(CHUNK)
    data_chunk=array('h',data)
    """
    Filter your signal!
    Afterwards, it doesn't contain lower frequencies anymore
    """
    signal_filtered, filter_state = signal.lfilter(b, a, data_chunk, zi=filter_state)

    """
    Convert amplitude to square amplitude, which is proportional to power,
    which is proportional to volume.
    """
    power = numpy.multiply(signal_filtered, signal_filtered)
    
    """
    Now, we'll also need to filter this – but this time with 
    a low-pass filter. We'll use the second filter we've designed:
    """
    power_filtered, filter_state = signal.lfilter(b_lpf, a_lpf, power, zi=filter_state_lpf)

    """
    Fun fact! We've now reduced the *actually occupied spectrum*
    significantly at the output of this rather hefty filter.
    That means of the 22.050 kHz that this signal originally could
    represent (says Nyquist), only maybe 1/200 actually still 
    contain 'change'. I.e. we can now more than safely simply ignore
    99 out of 100 samples - the others are definitely just nice and 
    smoothly interpolated between these.
    Now. That saves *a lot* of calculation. We only need to find the
    maximum on that subset! The real maximum power can't be much
    higher than what we see there.
    """
    maximum_power = max(power_filtered[::100])
    if maximum_power >= THRESHOLD:
        left_to_record = RECORD_SECONDS * RATE // CHUNK * CHUNK
        # I'm **REALLY** not versed with wavfile, and
        # I think you shouldn't be using that, but pysndfile
        # instead to write losslessly compressed (FLAC) audio
        # instead. Anyways, I haven't tried this, so debugging
        # is on you, sorry :)
        if not wavfile:
          timestring = time.strftime("%Y%m%d-%H%M%S")
          filename = f"RECORDING-{timestring}.wav"
          wavfile=wave.open(filename,'wb')
          wavfile.setnchannels(CHANNELS)
          wavfile.setsampwidth(audio.get_sample_size(FORMAT))
          wavfile.setframerate(RATE)
    
    if left_to_record > 0:
        """
        Write available full chunks!
        the // CHUNK * CHUNK trick just rounds down to the
        next multiple of CHUNK
        """
        
        write_now = (len(leftover_data) + len(data)) // CHUNK * CHUNK
        write_now = min(write_now, left_to_record)
        write_from_leftover = min(write_now, len(leftover_data))
        if write_from_leftover:
          wavfile.writeframes(bytes(leftover_data[:write_from_leftover]))
          leftover_data = [write_from_leftover]
          left_to_record -= write_from_leftover
          write_now -= write_from_leftover
        wavfile.writeframes(bytes(data[:write_now]))
        left_to_record -= write_now
        if left_to_record < CHUNK:
          wavfile.close()
          wavfile = None
        else:
          leftover_data = data[:write_now]

Obvious bugs:

The highpass and the low pass filters have both what we call group delay. So, the trigger comes a bit after the sounds have started. Is that bad? No, I think that would be less than 50 ms in delay, which you're probably OK with. And also, we don't actually capture 5s at the trigger, but 5s at the start of the frame, so on average we'll actually be a bit early.

This is, indeed, intended to "keep the recording active" for another 5s if a new super-threshold volume occurs within one recording time. I'm not sure that's your intention, but I think it is.


¹ Design parameters:

  • Highpass
  • IIR (because we don't care about phase, or niceness, this is just for detecting activity, so it should rather be low-effort! IIRs need less computation for the same amount of steepness, but they are not as nice as FIRs to the signal)
  • Elliptic (Good tradeoff between ability to quickly suppress a stopband, and being able to still have a reliable passband)
  • minimum order (we want the "cheapest" filter we get with these specifications, not the best filter we get for the "price" we pay in computational complexity)
  • Freq Unit kHz, $f_s=44.1$,
  • $f_{SB}$ (frequency where the stop band ends) 8
  • $A_{SB}$ the least attenuation for the stop band (in decibel): 50
  • $f_{PB}$ (frequency where the pass band starts) 11
  • $A_{PB}$ the ripple we allow in the passband, i.e. how un-flat the passband is allowed to be (decibel): 2.85
  • $f_C$ 11

² the second filter, the lowpass, looks like this:

low-pass filter magnitude response

$\endgroup$
6
  • $\begingroup$ Wow, thanks. I'll need some time to play around with this. I'll let you know if it works $\endgroup$
    – Thomas
    Jul 14 at 21:10
  • $\begingroup$ Just got back up, wandered to the PC and fixed another bug that occurred to me as I was pondering the night away. $\endgroup$ Jul 14 at 22:39
  • $\begingroup$ I've got as far as this line signal_filtered, filter_state = signal.lfilter(b, a, zi=filter_state) but get an error: TypeError: lfilter() missing 1 required positional argument: 'x' $\endgroup$
    – Thomas
    Jul 15 at 14:12
  • $\begingroup$ oops! That should be signal.lfilter(b, a, x=data_chunk, zi=filter_state)! $\endgroup$ Jul 15 at 15:46
  • $\begingroup$ Fixed that in the answer, too. $\endgroup$ Jul 15 at 15:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.