I've taken the following link http://coding-geek.com/how-shazam-works/#comments as a project to start practising and learning signal processing and python programming for me to get started in speech recognition and so on.

My experience with signal processing and python is no more then one week of heavy research on multiple papers, Wikipedia links and responses on stack and so on... thus there is not so much to expect from me and if there are issues please correct me, I would greatly appreciate.

I filled the whole code with comments as how I understand the whole process, that likely will help you somehow to get through it.

As I followed the article, I've started with initialising the parameters, all the songs(studio version, not live versions, thus they don't have noise) that each has a sampling rate 44.1khz and are in stereo mode:

from scipy.io import wavfile
from scipy.signal import butter, lfilter

import pyaudio
import math
import wave
import numpy as np
import matplotlib.pyplot as plt
import os

chunk = 1024

songs = [

notebook_path = os.path.abspath("Sound Recognition.ipynb")
file_path = os.path.join(os.path.dirname(notebook_path), "songs\\test\\" + songs[2] + ".wav")

# open up a wave
wf = wave.open(file_path, 'rb')

# sample width and framerate and number of channels
sample_width = wf.getsampwidth()
sampling_rate = wf.getframerate()
number_frames = wf.getnframes()
sample_channels = wf.getnchannels()

sampling_freq = 0
bin_freq = 0

# hamming window 1024 samples
window = np.hamming(chunk)

def sampling_details(chunk, sampling_rate):
    global sampling_freq
    sampling_freq = (chunk / sampling_rate)
    print("Sampling at a frequency of: %.4f hz" % sampling_freq)

    global bin_freq
    bin_freq = (sampling_rate / chunk)
    print("Frequency resolution for each frequency bin: %.4f hz" % bin_freq)

sampling_details(chunk, sampling_rate)

Here, I proceed to reading from audio source and parsing the data, after unpacking I proceed to converting the audio source from stereo to mono, and I'm starting with applying a low-pass filter where I cut frequency that exceeds the 1 / 4 of original sampling rate - 11 025hz and proceed to down sampling the audio source by taking the average of 4 samples that makes up each chunk as said in the article.

Then, I am applying the window function which I would expect to solve the problem of having too powerful upper bins and then getting the DFT result by applying FFT algorithm to each chunk and extracting the magnitude given data of complex numbers.

# opening stream
p = pyaudio.PyAudio()
stream = p.open(format = p.get_format_from_width(sample_width),
                channels = sample_channels,
                rate = sampling_rate,
                output = True)

# storing the data
out_data = []

# multichannel conversion to mono
def multi_to_mono(data):
    output = np.zeros(int(len(data) /  sample_channels))
    for c in range(sample_channels):
        # each nth is n channel where data is [0 1 .. n 0 1 .. n 0 1 .. n ...] 
        output += list(data[c: :sample_channels])

    output /= sample_channels

    return output.tolist()

# DOWNSAMPLING PART ***************************

# downsampling factor
downsampling_factor = 4

# new sampling rate after downsampling
sampling_rate_new = sampling_rate / downsampling_factor

# new hamming window with 512 samples
window = np.hamming(chunk / downsampling_factor)

# new sample len
new_chunk = int(chunk / downsampling_factor)

sampling_details(new_chunk, sampling_rate_new)

def downsampling(data):
    return (data[: 256] + data[256 : 512] \
        + data[512 : 768] + data[768 : 1024]) / 4

# FILTERING PART ***************************

# order of the filter
order = 6

nyquist_freq = 0.5 * sampling_rate

normal_cutoff = sampling_rate_new / nyquist_freq

b, a = butter(order, normal_cutoff, btype = 'low', analog = False)

# low pass filter
def low_pass_filter(data):
    return lfilter(b, a, data)

# half output
S = int(new_chunk * sample_width / 2)

# amplitude extractor from the dft data by taking the magnitude of a complex number 
magnitude_extractor = lambda x : np.sqrt(np.real(x) ** 2 + np.imag(x) ** 2) / S

# READING PART ***************************

# reading data
data = wf.readframes(chunk)

# playing stream and applying fft for each chunk
while len(data) == chunk * sample_width * 2:
    # write data out to the audio stream

    # unpack the data and parsing it
    unpacked_data = wave.struct.unpack("%dh" % (len(data) / sample_width), data)

    # conversion from multichannel input to mono output
    mono_data = multi_to_mono(unpacked_data)

    # filtering and downsampling:
    # filtering the data
    filtered_data = low_pass_filter(mono_data)

    # downsampling the data
    downsampled_data = downsampling(filtered_data)

    # times by the hamming window
    in_data = downsampled_data * window

    # taking the fft of the input signal, taking half as the latter is uncessary
    fft_output = np.fft.fft(in_data)[:S]

    # extracting the amplitude
    fft_data = magnitude_extractor(fft_output)

    # appending gathered data
    out_data += fft_data.tolist()

    # read more data that remains
    data = wf.readframes(chunk)

if data:


Then, I generate the logarithmic bands to filter the frequency bins to find the most potent ones:

Band 0 - from 0.000 hz to 73.403 hz

Band 1 - from 73.403 hz to 150.198 hz

Band 2 - from 150.198 hz to 307.339 hz


Band 7 - from 5387.985 hz to 11025.000 hz

Thus, each chunk that I sampled from audio source, I multiply by the frequency resolution that I've calculated way up that should remain the same even after applying a low pass filter and taking it through down sampling process, and taking the maximum ones that I use to calculate the average that will become a threshold that will tell what frequency band is the most prominent.

# offset_band because logarithm derivative is low in the beginning
offset_band = 6
count_band = 8
ceil_band = (offset_band + count_band)

# we substract one to get even 8 band ranges
base = np.power(sampling_rate_new, 1 / (ceil_band - 1))

# store band' frequency coverage
bands = [ base ** x for x in range(offset_band, ceil_band - 1) ]

message = "Band %d - from %.3f hz to %.3f hz"

print(message % (0, 0, bands[0]))
[ print(message % ((i - offset_band + 1), base ** i, base ** (i + 1))) for i in range(offset_band, (ceil_band -1)) ]

# getting the corresponding band based on freq
min_bin = base ** offset_band

def get_band(value):
    if value > min_bin:
        return math.floor(math.log(value, base) - offset_band + 1)
    else: return 0

class Band:
    def __init__(self, index, mag):
        self.index = index
        self.mag = mag

# we compute the strongest bins for the bands
output_bands = []

for i in range(0, len(out_data), new_chunk):
    max_band = np.zeros(count_band)

    for g in range(new_chunk):
        current_band = get_band(g * bin_freq)

        max_band[current_band] = max(max_band[current_band], out_data[g + i])

    # compute the average for the list of bands
    average = np.average(max_band).tolist()

    # filter the bands the exceed the average value
    output_bands_filtered = [ Band(index, mag) for index, mag in enumerate(max_band) if mag > average ]


output_freq = [ [ base ** (band.index + offset_band)  for band in bands ] for bands in output_bands]

down_scaler = 50
for index, values in enumerate(output_freq[: : down_scaler]):
    plt.scatter([index * sampling_freq * down_scaler] * len(values), values)

plt.xlabel("Time in s")
plt.ylabel("Frequency - hz")

And here is the filtered spectrogram, but as you see the most powerful bands are the ones at the start and at the end, which shouldn't be possible, I've used various window functions like hamming, blackman and (kazier with different parameters) that if you fit a polynomial to it you'll have low derivatives for very low values from start and end with slow increase of it till middle but the values taken from FFT algorithm are always too high no matter what.

I did some debugging and even if the values are very low at start and end after applying the window, the FFT algorithm returns extremely high values for complex numbers, both real and the imaginary part.

enter image description here

What I am doing wrong?


Your Answer

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

Browse other questions tagged or ask your own question.