Apply FIR filters to Extract Stego from WAV data

Question

I appreciate anyone that takes a moment to help me with this problem. I've been banging my head against the keyboard for a while, searching forums and DSP tutorials and I can't figure this problem out.

Here's a couple links I found very helpful from this forum on the subject: http://www.dspguide.com/ http://wiki.scipy.org/Cookbook/ApplyFIRFilter

So I have a WAV file with a 28 sec song/sound that I am trying to convert to frequency range, apply an FIR filter to isolate the 14.5 - 16kHz range. The WAV has a sampling rate of 32 kHz and contains an intermittent buzz that becomes clearer when I drop out all other frequencies besides the 14.5 -16kHz bands.

I have written this sample of python code, and believe I am close to isolating the signal. Its hard to tell.

# read data from file name fn, using low/high for bandpass
# low and high frequency cutoffs, and a 
def filter_frequencies(fn,low=0,high=0):
    data = wave.open(fn,'rb'); rate = data.getframerate()
    ntap = (rate / 4) + 1; nyq = rate / 2
    low,high = low/nyq,high/nyq
    chunks = int(data.getnframes()/(rate*4))
    ns = b''
    if low and high == 0: # lowpass
        bp = low; pz = True
    elif high and low == 0: # highpass
        bp = high; pz = False
    else: # bandpass
        bp = [low,high]; pz = False

    b = firwin(ntap, bp, pass_zero=pz, window='blackman', nyq=nyq)
    y = b''
    for num in range(chunks):
        # read 4 secs of frames from wav file
        seg = np.fromstring(data.readframes(rate*4),dtype=np.int8)
        #seg = np.fromstring(data.readframes(data.getnframes()),dtype=np.int8)
        # convert frame amplitudes to their frequency values
        freq = np.fft.rfft(seg)
        # perform lowpass, highpass or bandpass filter
        cr = fftconvolve(freq,b,mode='valid')
        # converts frequencies to binary string
        y += cr.tostring()
    data.close()
    return y

I am a noob to DSP, so again I would really appreciate your help on if I am doing the FFT to inverse FFT correctly.

Below is the code I am using to transform the output from the above function back into WAV bytes format.

def frequency_to_bytes(freq):
    # invert the Fast Fourier Transform
    nf = np.fft.irfft(freq)
    # convert frequencies back to integer values
    bs = np.ravel(nf).astype(np.int8)
    return bs.tostring()

P.s. if you have seen this problem before, or know what challenge I am trying to solve, please no spoilers. I have been working on this for a week solid.

edit: added some background links I have been using as reference.

Answer 1

I never worked with Python or SciPy, but I don't see the sense in the actions you apparently perform in your code. firwin provides you with the time-domain filter coefficients, which are to be convolved with the time-domain signal you wish to filter. What you do in your code, is convolve the time-domain filter coefficients with the frequency-domain signal - which makes absolutely no sense.

What is the point of converting the signal to freq-domain (doing FFT)? Is the stego analysis/decoding performed in time- or frequency-domain?

If your processing chain requires conversion of the signal to freq-domain, for some purpose, then you don't really need the FIR filter, since you can get rid of the undesired frequency elements once you have the signal in the frequency domain, by simply zeroing-out the undesired frequency bins of the signal.

If the only reason for calculating FFT of the signal is to filter it - then there is no need to do it, as you can simply convolve the time-domain filter coefficients with the time-domain signal, as I explained.

Answer 2

The amazing engineering field called Digital Signal Processing (DSP) comes from a simple analysis of discrete-time signals to complex adaptive filters. But, to start, a nice idea is to think of the discrete-time signals as a vector, where each element of this vector is a sampled value of the original, continuous-time signal. Once you get the samples in a vector form, you can apply different digital signal techniques to this vector.

Unfortunately, on Python, moving from audio files to NumPy array vector is rather cumbersome, as you could notice... If you don't idolize one programming language over other, I highly suggest trying out MatLab/Octave. Matlab makes the samples access from files straightforward. audioread() makes this task to you :) And there are a lot of toolboxes designed specifically for DSP.

Nevertheless, if you really intend to get into Python for this, I'll give you a step-by-step to guide you.

---

1. Get the samples

The easiest way the get the samples from the .wav file is:

from scipy.io import wavfile

sampling_rate, samples = wavfile.read(f'/path/to/file.wav')

---

Alternatively, you could use the wave and struct package to get the samples:

import numpy as np
import wave, struct

wav_file = wave.open(f'/path/to/file.wav', 'rb')
# from .wav file to binary data in hexadecimal
binary_data = wav_file.readframes(wav_file.getnframes())
# from binary file to samples
s = np.array(struct.unpack('{n}h'.format(n=wav_file.getnframes()*wav_file.getnchannels()), binary_data))

Answering your question: binary_data is a bytes object, which is not human-readable and can only make sense to a machine. You can validate this statement typing type(binary_data). If you really want to understand a little bit more about this bunch of odd characters, click here.

If your audio is stereo (that is, has 2 channels), you can reshape this signal to achieve the same format obtained with scipy.io

s_like_scipy = s.reshape(-1, wav_file.getnchannels())

Each column is a chanell. In either way, the samples obtained from the .wav file can be used to plot and understand the temporal behavior of the signal.

In both alternatives, the samples obtained from the files are represented in the Linear Pulse Code Modulation (LPCM)

---

2. Do digital signal processing stuffs onto the audio samples

I'll leave that part up to you :) But this is a nice book to take you through DSP. Unfortunately, I don't know good books with Python, they are usually horrible books... But do not worry about it, the theory can be applied in the very same way using any programming language, as long as you domain that language.

Whatever the book you pick up, stick with the classical authors, such as Proakis, Oppenheim, and so on... Do not care about the language programming they use. For a more practical guide of DPS for audio using Python, see this page.

3. Play the filtered audio samples

import pyaudio

p = pyaudio.PyAudio()
stream = p.open(format = p.get_format_from_width(wav_file.getsampwidth()),
                channels = wav_file.getnchannels(),
                rate = wav_file.getframerate(),
                output = True)
# from samples to the new binary file
new_binary_data = struct.pack('{}h'.format(len(s)), *s)
stream.write(new_binary_data)

where wav_file.getsampwidth() is the number of bytes per sample, and wav_file.getframerate() is the sampling rate. Just use the same parameters of the input audio.

---

4. Save the result in a new .wav file

wav_file=wave.open('/phat/to/new_file.wav', 'w')

wav_file.setparams((nchannels, sampwidth, sampling_rate, nframes, "NONE", "not compressed"))

for sample in s:
   wav_file.writeframes(struct.pack('h', int(sample)))

where nchannels is the number of channels, sampwidth is the number of bytes per samples, sampling_rate is the sampling rate, nframes is the total number of samples.