Quantum Computing to find frequencies in an audio file

Because why not...

Do read this: article by BCG predicting the economic impact of quantum computing in the NISQ era and beyond - which estimates the economic impact of quantum computing in the next 5 years to be $2 - 5 billion.

My inner quantum physicist sprung into action and I decided to make a tangible application on a quantum computer - however simple it maybe.

All code can be found at this repository

So What’s Happening Exactly ?

First things first - what am I trying here ? - basically I am trying to build a frequency detector on a quantum computer, audio file goes in and most dominant frequencies come out, thats it. (some might ask - aren’t you making a glorified version of a quantum version of discrete fourier transform ? -> well yeh but sshhh!). I plan to start with simple audio signal of a pure sine wave. If things work out I may try to record a note played on a guitar.

TLDR; - I want to make a frequency detector using a quantum computer

Why should I care ?

Apart from having the potential to crack your facebook/instagram/bank account passwords, QFT or Fourier transform is one of the central tools in mathematics and is used almost in every piece of tech that you see around you in some way. The difference between Quantum Fourier Transform and classical Fast Fourier transform is in the speed and in the way the data is represented physically. Classically, a dimension vector would need floating point numbers. On a quantum computer, the QFT operates on the wave function needing only qubits, exponentially saving space. The best classical FFT runs in time and the QFT runs in time where again is the dimension of the vector. Also the classical FFT must take time to even read the input!

The vanilla QFT algorithm takes quantum gates, but there are very efficient approximate versions which need only gates possibly giving more speedup.

1. Quantum Fourier Transform

QFT is the quantum version of the discrete fourier transform. In QFT we encode the input vector as a set of amplitudes of the basis states of the quantum system. Hence, the number of points we can operate over is restricted by the number of qubits available to us.

The classical fourier transform acts on a vector and calculates where

where

The quantum fourier transform does something similar. It acts on a quantum state

and maps it to a quantum state

according to the equation:

I am using the convention which makes the QFT have the same effect as the inverse discrete fourier transform. Depending on the implementation of the circuit this convention can change.

If we deal with basis states then QFT can also be expressed as a map:

If you are the sort of person who better understands stuff in terms of unitaries, the equivalent unitary for QFT is

QFT Unitary

1.1 Building the quantum circuit for QFT

* Quirk visualizer by Craig Gidney

The quantum processors that we currently have encode operations in quantum gates. The next few lines rewrite the QFT equation into a form which can be implemented as a quantum circuit using the commonly available quantum gates. You can skip this if you want.

1.2 Digression

To understand the above equation it would help to expand it using a small n. For example for a 2 qubit system, the basis states are

1.3 Moving on

Hence,

Hence, the original QFT equation can be rewritten in terms of individual qubits as :

Now,

The exponent can be further simplified as:

where are the binary components of

Notice that above is always a whole number, and

where the fractional binary notation as below is used:

Hence,

Finally we can write:

Now that we have the QFT formula as a tensor product of each individual qubit, it would be easier to see what gates would get us the superposition in the above equation. The QFT circuit mainly makes use of two gates - Hadamard gate (H) and a parametric controlled phase gate .

  1. Hadamard Gate (H)

Hadamard Gate

The Hadamard gate is the most commonly used gate to put a pure qubit into a superposition state. It coverts qubits from the Pauli-Z basis to the Pauli-X basis. Its action is given as follows:

and in general:

  1. Controlled Phase Gate

Controlled phase gate

This is a two qubit controlled gate assuming that the control qubit is the most significant qubit, the action of this gate can be expressed as:

Using these two gates, and the equation derived for QFT, we can construct the full circuit to get us the desired superposition like so

* Wikipedia

There’s a caveat here, the above circuit reverses the qubit order upside down, so many implementations swap the qubits to restore the qubit order. Depending on the application this swap may or may not be necessary. If you are having trouble understanding how does the circuit give us the needed superposition, I urge you to try writing it down for the case of 2 qubits. That’s what I did.

QFT circuit can be very easily made using Qiskit :

def qft_rotations(circuit, n):
    """Performs qft on the first n qubits in circuit (without swaps)"""
    if n == 0:
        return circuit
    n -= 1
    circuit.h(n)
    for qubit in range(n):
        circuit.cu1(pi/2**(n-qubit), qubit, n)
    #note the recursion
    qft_rotations(circuit, n)

def swap_registers(circuit, n):
    for qubit in range(n//2):
        circuit.swap(qubit, n-qubit-1)
    return circuit

def qft(circuit, n):
    """QFT on the first n qubits in circuit"""
    qft_rotations(circuit, n)
    swap_registers(circuit, n)
    return circuit

2. Representing Audio in quantum data

An audio file is simply a bunch of samples sampled at a specific sampling rate. Usually this rate is 44100Hz for CD quality audio. This means that the analog waveform is measured at 44100 equally spaced points within a second.

900Hz sine wave
900Hz sine wave

If we want to use QFT on such a wave we need to encode the value of each sample into the amplitudes of the basis states of our system. Consider this 3-qubit example:

3-qubit samples
3-qubit samples

For a 3 qubit system, we can afford samples, which are denoted by the Cs in the above figure,

Note that Since is a quantum state, we need

Therefore, the Cs we sample need to be normalised before we can make a quantum state out of them.

3. Arbitrary state preparation

After getting all the Cs, its not trivial to put a system into the exact superposition that we want. To do that, you want to know what Schmidt decomposition is

3.1 Schmidt Decomposition Theorem

Any vector can be expressed in the form

for non-negative real and orthonormal sets and . There are density operators on and on such that :

for all observables and on and respectively, and the can be chosen to be the eigenvectors of corresponding to the non-zero eigenvalues , the vectors , the corresponding eigenvectors for , and the positive scalars

3.2 State preparation

A general n-qubit state is fully described by real parameters. These parameters are introduced sequentially in the form of C-NOT gates and single qubit rotation unitaries U. For example ,we can look at the process for a 4 qubit system below. For the general case I’d refer you to the paper “Quantum State Preparation with Universal Gate Decompositions”

In the case of four qubits, the Hilbert space can be factorized into two parts each of two qubits. An arbitrary pure state of four qubits can then be represented using the standard Schmidt decomposition as

After that the entire process can be broken down into four phases.

Phase 1

Starting with the initial state of all zeros, we generate the state with generalized schmidt coefficients on the first two qubits:

This can be done using a single C-NOT gate and single qubit rotations. For more details about this step refer “Decompositions of general quantum gates”

Phase 2

We perform two C-NOT operations, one with the control on the first qubit and the target on the third qubit and the other one with the control on the second qubit and the target on the fourth qubit. In such a way we can “copy” the basis states of the first two qubits onto the respective states of the second two qubits. In this way we obtain a state of four qubits, which has the same Schmidt decomposition coefficients as the target

Phase 3

Keeping the Schmidt decomposition form we apply the unitary operation that transforms the basis states of the first two qubits into the four states, we obtain:

Phase 4

In the final phase of the circuit we perform a unitary operation on the third and fourth qubit in order to transform their computational basis states into the Schmidt basis states

This completes the state preparation.

3.3 State preparation in Qiskit

Thankfully, qiskit provides us with inbuilt functions for preparing arbitrary quantum states so we don’t have to do the above mathematical fuckery ourselves. Using Qiskit, a state preparer can be easily implemented using the following snippet of code:

def prepare_circuit(samples, normalize=True):

    """
    Args:
    amplitudes: List - A list of amplitudes with length equal to power of 2
    normalize: Bool - Optional flag to control normalization of samples, True by default
    Returns:
    circuit: QuantumCircuit - a quantum circuit initialized to the state given by amplitudes
    """
    num_amplitudes = len(samples)
    assert isPow2(num_amplitudes), 'len(amplitudes) should be power of 2'

    num_qubits = int(getlog2(num_amplitudes))
    q = QuantumRegister(num_qubits)
    qc = QuantumCircuit(q)

    if(normalize):
        ampls = samples / np.linalg.norm(samples)
    else:
        ampls = samples

    qc.initialize(ampls, [q[i] for i in range(num_qubits)])

    return qc

This is an example of how it looks, notice that it only uses C-NOT gates and single qubit rotations

Init Circuit

3.4 Important note

Given the noisy quantum processors that we currently have, its imperative to minimise the number of gates we have as much as possible. In my experience, in a 4-qubit IBM-Q system, as few as 50 gates was enough to introduce noise enough to get garbage output.

Hence, practically you might want to do your own optimizations to the state preparation circuit (e.g. getting rid of small angle rotations), which the qiskit inbuilt function does not allow.

4. Running the experiment on a quantum simulator

Qiskit provides many simulator backend options, here I am going to use the qasm simulator backend, which mimics the behaviour of a noise-free real backend.

Since its a simulator, I am working with 16 qubits, which means I will be able to load in samples of my audio file.

The entire circuit can be split into 4 parts -

  1. WAV file processor
  2. Quantum state preparation
  3. Quantum Fourier Transform
  4. Measurement

I have used the pydub package to process the wav file. All quantum operations are done using Qiskit. The full code can be found in the linked github repository.

Pipeline

Step 1 - Input Audio file and prepare a quantum state

n_qubits = 16
audio = AudioSegment.from_file('900hz.wav')
samples = audio.get_array_of_samples()[:2**n_qubits]
plt.xlabel('sample_num')
plt.ylabel('value')
plt.plot(list(range(len(samples))), samples)

# prepare circuit
qcs = prepare_circuit_from_samples(samples, n_qubits)

Step 2 - Apply QFT and get amplitudes in the fourier space

qcs = prepare_circuit_from_samples(samples, n_qubits)
qft(qcs, n_qubits)
qcs.measure_all()

qasm_backend = Aer.get_backend('qasm_simulator')
out = execute(qcs, qasm_backend, shots=8192).result()

Step 3 - Measurement

When a quantum circuit is run, at the end each qubit is measured in the Pauli-Z basis. Hence, all the qubits collapse into one of the two eigenstates of the Pauli-Z operator. But they do so with the probabilities dictated by the superposition that we put the qubits in. The big mathematical fuckery that we did in section 1 was precisely to get each qubit in a superposition that would make the qubits collapse in a way that gives us the fourier transform of the input amplitudes.

If the final output of the state is given by :

where are the basis states, then the information we need is in the

However since we can only observe the distribution of the final qubits, we can only measure

To get this distribution we run the circuit many times, which is controlled by the shots parameter.

After measurement, we still need to convert the distribution that we get into actual frequencies. This is a pretty standard operation on discrete fourier transforms. If you feel uncomfortable understanding the next few lines I’d suggest you see the discrete fourier transform page on wikipedia.

out = execute(qcs, qasm_backend, shots=8192).result()
counts = out.get_counts()
fft = get_fft_from_counts(counts, n_qubits)[:n_samples//2]
plot_samples(fft[:2000])

f-space

We see a nice peak at about 1300. To get the value of frequency we have to multiply this peak by . Note that we can only measure upto - which is the Nyquist limit. If you don’t know what that means google is your friend.

top_indices = np.argsort(-np.array(fft))
freqs = top_indices*frame_rate/n_samples
# get top 5 detected frequencies
freqs[:5]
=> Prints out: array([899.68414307, 900.35705566, 899.01123047, 901.02996826,
       898.33831787])

Nice! Great Success

Running the experiment on real quantum hardware!

Totally different beast that is. Real quantum hardware that we currently have is extremely noisy.

The way you access real quantum hardware is by making use of IBM - Q Experience. Currently it offers various real quantum backends. One of the backends has a capacity of 16 qubits (ibmq_16_melbourne), while the rest have only 5 qubits

Trial 1

I tried the circuit that I built in the previous sections as is on the 16 qubit backend - it failed.

No surprises there, I was sure it couldn’t be that easy. The problem was that I was exceeding the max depth (number of gates) allowed by the backend.

error1

Result: Failure

Trial 2

I figured perhaps the 16 qubit experiment is a bit too much, so I tried simplifying it to 5 qubits, with a very low sample size of samples. With such a small number of samples, to get a decent frequency resolution, I had to drastically decrease the sampling rate to , which gave a frequency resolution of

I first ran the experiment on the qasm_simulator to see how the desired output should look like:

n_qubits = 5
n_samples = 2**n_qubits
audio = AudioSegment.from_file('900hz.wav')
audio = audio.set_frame_rate(2000)
frame_rate = audio.frame_rate #2000
samples = audio.get_array_of_samples()[:n_samples]
plt.xlabel('sample_num')
plt.ylabel('value')
plt.plot(list(range(len(samples))), samples)

sine_900_2000

With so few samples, the plot starts loosing its sine-wavy looks.

qcs = prepare_circuit_from_samples(samples, n_qubits)
qft(qcs, n_qubits)
qcs.measure_all()

qasm_backend = Aer.get_backend('qasm_simulator')
out = execute(qcs, qasm_backend, shots=8192).result()
counts = out.get_counts()
fft = get_fft_from_counts(counts, n_qubits)[:n_samples//2]


plt.xlabel('sample_num')
plt.ylabel('value')
plt.plot(list(range(len(fft[:]))), fft[:])

top_indices = np.argsort(-np.array(fft))
freqs = top_indices*frame_rate/n_samples
# get top 5 detected frequencies
freqs[:5]
=> Prints: array([875. , 937.5, 812.5, 750. , 687.5])

f32

The plot makes sense with one strong peak. Alright lets get real and try it on real hardware.

provider = IBMQ.get_provider(hub='ibm-q')
backend = least_busy(provider.backends(filters=lambda x: x.configuration().n_qubits >= n_qubits
                                       and not x.configuration().simulator
                                       and x.status().operational==True))
print("least busy backend: ", backend)
=> least busy backend:  ibmqx2
out = execute(qcs, qasm_backend, shots=8192).result()

You can see the results of your experiments in the IBM Q Experience portal. For this particular experiment, it showed me this:

ibmq_32

Which is gibberish, we ideally should get one single peak. So something’s clearly wrong.

Result: Failure

Trial 3

I thought maybe there is something wrong in the way the state is prepared from the audio file. This was supported by the fact that more than half of the gates in the circuit were actually used for state preparation.

I reduced the number of qubits to 4, and created a dummy state comprising of basis states with linearly increasing amplitudes.

n_qubits = 4
dummy_state = [i for i in range(2**n_qubits)]
=> dummy_state = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]

qcs = prepare_circuit_from_samples(dummy_state, n_qubits, normalize=True)
qcs.measure_all()

The circuit as visualized on IBM Q portal:

After running on a simulator, the results look like :

f_16

Which looks good as expected.

After running it on real hardware:

ibmq_16_dummy_histogram.png

Its noisy again! borderline gibberish. Ideally the amplitudes should go up linearly. But at least the problem has been isolated.

Lets try the same dummy state preparation with 3 qubits.

n_qubits = 3
dummy_state = [i for i in range(2**n_qubits)]
=> dummy_state = [0, 1, 2, 3, 4, 5, 6, 7]

qcs = prepare_circuit_from_samples(dummy_state, n_qubits, normalize=True)
qcs.measure_all()

The circuit as visualized on IBM Q portal:

After running on a simulator, the results look like :

f_16

Which looks good as expected.

After running it on real hardware:

ibmq_8_dummy_histogram.png

Now we are talking! This looks much more acceptable than previous experiments. Which sadly means that we will have to restrict ourselves a a meagre 3 qubits, which is only 8 samples. Lets see what we can do with it.

Result: almost Success

Trial 4

Running the entire circuit, with a sampling rate of , and only 3 qubits

n_qubits = 3
n_samples = 2**n_qubits
audio = AudioSegment.from_file('900hz.wav')
audio = audio.set_frame_rate(2000)
frame_rate = audio.frame_rate #2000
samples = audio.get_array_of_samples()[:n_samples]

qcs = prepare_circuit_from_samples(samples, n_qubits)
qft(qcs, n_qubits)
qcs.measure_all()

provider = IBMQ.get_provider(hub='ibm-q')
backend = least_busy(provider.backends(filters=lambda x: x.configuration().n_qubits >= n_qubits
                                       and not x.configuration().simulator
                                       and x.status().operational==True))
print("least busy backend: ", backend)
=> least busy backend:  ibmqx2
out = execute(qcs, qasm_backend, shots=8192).result()

The result looks like:

ibmq_8_histogram.png

Nice! we see a clear peak. The peak seems to be at which is . To get the associated frequency :

Nice!

Result: Success! (?)

I gotta say..

This work is more like a proof of concept. It may not look like much for now but as we move towards fault tolerant, less error prone quantum hardware, our capacities would increase exponentially.

There could be a lot of ways in which this work can be improved even with the current hardware. For starters:

  1. No error correction was used.
  2. Approximate quantum fourier transform could’ve been used, where we get rid of small angle rotations to minimize the number of gates.
  3. State preparation routine could be made approximate to further minimize number of gates.

References:

  1. Code for this project
  2. BCG Article talking about impact of quantum computers
  3. What is NISQ Era quantum hardware?
  4. “Quantum State Preparation with Universal Gate Decompositions”
  5. “Decompositions of general quantum gates”
  6. “Approximate Quantum Fourier Transform”
  7. “Quirk quantum circuit visualizer”