Abstract: Fourier Transform is a powerful concept, used in a variety of fields, from pure mathematics to audio engineering and even finance.

Scipy. FFT module

Fourier transform is an important tool in many applications, especially in scientific computing and data science. As a result, SciPy has long provided its implementation and associated transformations. Initially, SciPy provided the SciPy. Fftpack module, but they later updated their implementation and moved it into the SciPy. FFT module.

SciPy is packed with features. For a more general introduction to the library, see Scientific Python: Optimization with SciPy.

Install SciPy and Matplotlib

To begin, you need to install SciPy and Matplotlib. You can do this in one of two ways:

  1. ** Install using Anaconda: ** Download and install Anaconda Individual Edition. It comes with SciPy and Matplotlib, so once you follow the steps in the installer, you’re done!
  2. If you have PIP installed, you can install the library using the following command:
$ python -m pip install -U scipy matplotlib
Copy the code

You can verify that the installation works by typing Python into the terminal and running the following code:

\>>> \>>> import scipy, Matplotlib \ > > > print (scipy. \ _ \ _file \ _ \ _)/usr/local/lib/python3.6 / dist - packages/scipy / \ _ \ _init \ _ \ _ py \ > > > Print (matplotlib. \ _ \ _file \ _ \ _)/usr/local/lib/python3.6 / dist - packages/matplotlib / \ _ \ _init \ _ \ _ pyCopy the code

This code imports SciPy and Matplotlib and prints the location of the module. Your computer may display a different path, but as long as it prints the path, the installation is successful.

SciPy is now installed! Now it’s time to look at the difference between scipy.fft and scipy.fftpack.

Scipy. FFT contrast scipy. Fftpack

When you look at the SciPy documentation, you might come across two modules that look very similar:

  1. scipy.fft

  2. scipy.fftpack

The scipy.fft module is relatively new and should take precedence over scipy.fftpack. You can read more about the changes in the release notes for SciPy 1.4.0, but here’s a quick summary:

  • Scipy.fft has an improved API.

  • Scipy. FFT allows the use of multiple workers, which can provide a speed boost in some cases.

  • Scipy. Fftpack is considered legacy and SCIpy recommends that Scipy. FFT be switched over.

Unless you have a good reason to use Scipy.fftpack, you should stick with scipy.fft.

Scipy. FFT contrast numpy. FFT

SciPy’s fast Fourier Transform (FFT) implementation contains more features and is more likely to fix errors than NumPy’s implementation. If you have a choice, you should use the SciPy implementation.

NumPy maintains the FFT implementation for backward compatibility, although the authors feel that features like Fourier transform are best placed in SciPy. See the SciPy FAQ for more details.

Fourier Transform

Fourier analysis is the study of the decomposition of mathematical functions into a series of simpler trigonometric functions. A tool in this field is the Fourier transform, which is used to decompose functions into their component frequencies.

Well, that’s a very dense definition. For the purposes of this tutorial, the Fourier transform is a tool that allows you to take a signal and see the power of each frequency within it. Take a look at the important terms in this sentence:

  • A signal is information that changes over time. For example, audio, video, and voltage traces are examples of signals.

  • A frequency is the rate at which something repeats itself. For example, a clock ticks at a frequency of one Hertz (Hz), or repeats once per second.

  • In this case, power represents only the intensity of each frequency.

Here is a visual illustration of the frequency and power of some sine waves:

The peaks of high frequency sine waves are closer together than those of low frequency sine waves because they repeat more frequently. The low power sine wave has a smaller peak than the other two sine waves.

To illustrate this more concretely, suppose you use the Fourier transform on a recording of someone playing three notes simultaneously on a piano. The resulting spectrum will show three peaks, one for each note. If the person plays one note more softly than the others, the intensity of that note’s frequency will be lower than the other two.

Here’s what the piano example looks like visually:

The highest note on the piano is quieter than the other two notes, so the spectrum of that note has a lower peak.

Why Fourier Transform?

The Fourier transform is useful in many applications. Shazam and other music-recognition services, for example, use The Fourier transform to identify songs.

JPEG compression uses a variant of the Fourier transform to remove the high frequency components of the image. Speech recognition uses Fourier transform and correlation transform to recover speech from the original audio.

In general, if you need to look at the frequencies in the signal, you need to perform a Fourier transform. If processing a signal in the time domain is difficult, it is worth trying to move it into the frequency domain using the Fourier transform. In the next section, you will learn about the difference between the time domain and the frequency domain.

Time Domain vs Frequency Domain

In the rest of this tutorial, you will see the terms Time domain and Frequency Domain. The terms refer to two different ways of looking at a signal, either its component frequency or information that changes over time.

In the time domain, the signal is a wave whose amplitude (Y-axis) varies with time (X-axis). You’re probably used to seeing graphs in the time domain, such as this one:

This is a picture of some audio, and it’s a time domain signal. The horizontal axis represents time and the vertical axis represents amplitude.

In the frequency domain, the signal is represented as a series of frequencies (the X-axis), each with an associated power (the Y-axis). The following figure shows the above audio signal after Fourier Transform:

Here, the previous audio signal is represented by its constituent frequency. Each frequency at the bottom has an associated power that produces the spectrum you see.

For more information about the frequency domain, see the DeepAI vocabulary entry.

The type of Fourier Transform

The Fourier transform can be subdivided into different types of transformations. The most basic refinement is based on the data type of the transformation operation: continuous function or discrete function. This tutorial will deal only with the discrete Fourier Transform (DFT).

Even in this tutorial, you will often see the terms DFT and FFT used interchangeably. However, they are not exactly the same. The ** fast Fourier Transform (FFT) ** is an algorithm used to compute the discrete Fourier Transform (DFT), and the DFT is the transform itself.

Another difference you’ll see in the scipy.fft library is the difference between different types of input. FFT () accepts complex numeric input, and RFFT () accepts real numeric input. Skip to using the Fast Fourier Transform (FFT) section to learn about complex numbers and real numbers.

Two other transformations are closely related to the DFT: the discrete cosine transform (DCT) and the discrete Sine transform (DST). You will learn about this in the discrete cosine and sine transforms section.

Practical example: Remove unwanted noise from audio

To help you understand the Fourier transform and what you can do with it, you’ll filter some audio. First, you’ll create an audio signal with a high-pitched hum, and then you’ll use the Fourier transform to remove the hum.

Create a signal

Sine waves are sometimes called pure tones because they represent a single frequency. You will use sine waves to generate audio because they will form different peaks in the generated spectrum.

Another advantage of sine waves is that they can be generated directly using NumPy. If you haven’t used NumPy before, what is NumPy?

Here is some code for generating sine waves:

import numpy as np  
from matplotlib import pyplot as plt  
  
SAMPLE\_RATE = 44100  # Hertz  
DURATION = 5  # Seconds  
  
def generate\_sine\_wave(freq, sample\_rate, duration):  
    x = np.linspace(0, duration, sample\_rate \* duration, endpoint=False)  
    frequencies = x \* freq  
    # 2pi because np.sin takes radians  
    y = np.sin((2 \* np.pi) \* frequencies)  
    return x, y  
  
\# Generate a 2 hertz sine wave that lasts for 5 seconds  
x, y = generate\_sine\_wave(2, SAMPLE\_RATE, DURATION)  
plt.plot(x, y)  
plt.show()
Copy the code

After you import with NumPy and Matplotlib, you can define two constants:

  1. SAMPLE_RATE determines how many data points a signal uses per second to represent a sine wave. Therefore, if the signal has a sampling rate of 10 Hz and is a sine wave of 5 seconds, it will have 10 * 5 = 50 data points.

  2. DURATION is the length of the sample to be generated.

Next, you define a function to generate the sine wave, because you will use it many times in the future. This function takes the frequency, freq and then returns the x and Y values used to plot the waveform.

The sine wave’s x coordinates are DURATION evenly distributed between 0 and, so the code uses NumPylinspace() to generate them. It requires a starting value, an ending value, and the number of samples to generate. Setting endpoint=False is important for the Fourier transform to work properly because it assumes that signals are periodic.

Np.sin () evaluates the sine function at each x-coordinate. The result is multiplied by the frequency so that the sine wave oscillates at that frequency, multiplying the product by 2π to convert the input value into radians.

** Note: ** If you haven’t studied much trigonometry before, or if you need a refresher, check out khan Academy’s trigonometry courses.

Once you define the function, you can use it to generate a 2-Hertz sine wave that lasts 5 seconds and draw it using Matplotlib. Your sine wave diagram should look like this:

The X-axis shows time in seconds, and since there are two peaks in each second, you can see the sine wave oscillating twice per second. The frequency of this sine wave is too low to hear, so in the next section, you will generate some higher frequency sine waves, and you will learn how to mix them.

Mixed audio signal

The good news is that mixing audio signals involves only two steps:

  1. Add up the signals

  2. Standardization result

Before mixing signals together, you need to generate them:

\_, nice\_tone = generate\_sine\_wave(400, SAMPLE\_RATE, DURATION)  
\_, noise\_tone = generate\_sine\_wave(4000, SAMPLE\_RATE, DURATION)  
noise\_tone = noise\_tone \* 0.3  
  
mixed\_tone = nice\_tone + noise\_tone
Copy the code

There is nothing new in this code example. It generates the mezzo and soprano noise_tones assigned to the variables nice_tone and respectively. You will use the treble as the noise you do not need, so it will be multiplied by 0.3 to reduce its power. The code then adds the tones together. Notice that you use the underscore (_) to discard the value generate_sine_wave() returned by x.

The next step is to normalize, or scale, the signal to fit the target format. Because of how you will store audio later, your target format is a 16-bit integer ranging from -32768 to 32767:

normalized\_tone = np.int16((mixed\_tone / mixed\_tone.max()) \* 32767)  
  
plt.plot(normalized\_tone\[:1000\])  
plt.show()
Copy the code

Here, the code scales mixed_tone to fit 16-bit integers, and then uses NumPy’s Np.int16. Divide mixed_tone at its maximum scale to -1 and between 1. When the signal is multiplied by 32767, it scales between -32767 and 32767, roughly in the range of np.int16. This code only draws the first 1000 sample so that you can see the structure of the signal more clearly.

Your plot should look something like this:

The signal looks like a distorted sine wave. The sine wave you see is the 400 Hz tone you generated, and the distortion is the 4000 Hz tone. If you look closely, you can see that the distortion takes the shape of a sine wave.

To listen to audio, you need to store it in a format that the audio player can read. The simplest way is to use SciPy’s wavFile.write method to store it in a WAV file. 16-bit integers are the standard data type for WAV files, so you need to normalize signals to 16-bit integers:

from scipy.io.wavfile import write  
  
\# Remember SAMPLE\_RATE = 44100 Hz is our playback rate  
write("mysinewave.wav", SAMPLE\_RATE, normalized\_tone)
Copy the code

This code will write to a file in mysinewave.wav in the directory where you ran the Python script. You can then listen to this file using any audio player or even Python. You will hear a lower pitch and a higher pitch. These are the 400 Hz and 4000 Hz sine waves that you mix.

When this is done, your audio sample is ready. The next step is to use the Fourier transform to remove the treble!

Using Fast Fourier Transform (FFT)

It’s time to use FFT on the generated audio. FFT is an algorithm that implements the Fourier transform and calculates the spectrum of a signal in the time domain, such as your audio:

from scipy.fft import fft, fftfreq \# Number of samples in normalized\_tone N = SAMPLE\_RATE \* DURATION yf = fft(normalized\_tone) xf = fftfreq(N,  1 / SAMPLE\_RATE) plt.plot(xf, np.abs(yf)) plt.show()Copy the code

This code will compute the Fourier Transform of the generated audio and draw it. Before breaking it down, take a look at the plot it spawned:

You can see the two peaks in the positive frequency and the mirror image of those peaks in the negative frequency. Positive frequency peaks are located at 400 Hz and 4000 Hz, which correspond to the frequency at which you input audio.

Fourier Transform has transformed your complex, weak signal into the frequency it contains. Because you only put in two frequencies, only two frequencies come out. Negative and positive symmetry is a side effect of putting real input into Fourier Transform, but you will hear more about this later.

In the first few lines, you import the function that scipy.fft will use later and define a variable, N, to store the total number of samples in the signal.

After this comes the most important part, computing the Fourier Transform:

yf = fft(normalized\_tone)  
xf = fftfreq(N, 1 / SAMPLE\_RATE)
Copy the code

The code calls two very important functions:

  1. FFT () evaluates the transformation itself.

  2. Fftfreq () computes the frequency FFT () of each bin center in the output. Without that, you can’t plot the X-axis on the spectrum.

A box is already grouped like a histogram in a range of values. For more information about bin, see this signal handling stack swapping issue. For the purposes of this tutorial, you can think of them as individual values.

Once you have the result values of the Fourier Transform and their corresponding frequencies, you can draw them:

plt.plot(xf, np.abs(yf))  
plt.show()
Copy the code

The interesting part of this code is what you do with YF before you draw. You call Np.abs () yF because its value is complex.

A complex number is a number with two parts, real and imaginary. There are many reasons why defining such numbers is useful, but all you need to know for now is that they exist.

Mathematicians usually write complex numbers in the form a + bi, where A is the real part and b is the imaginary part. The I after b means that b is imaginary.

** Note: ** Sometimes you’ll see complex numbers written with I, and sometimes you’ll see complex numbers written with j, such as 2 + 3i and 2 + 3j. It’s the same thing, but I is used more by mathematicians, and J is used more by engineers.

For more information on complex numbers, check out the Khan Academy courses or Math is Fun page.

Since the complex number has two parts, plotting them against frequencies on a two-dimensional axis requires you to compute a value from them. This is where Np.abs () comes in. It computes the square root of a squared plus b squared of the complex number, which is the total size of the two numbers, but the individual values are important.

** note: ** by the way with FFT (), you may have noticed that the maximum frequency returned is just over 20,000 hz, or 22050Hz, to be exact. This value, which is exactly half of our sampling rate, is called the Nyquist frequency.

This is a basic concept in signal processing, meaning that your sampling rate must be at least twice the highest frequency of the signal.

Make it faster RFFT ()

The spectrum of the output of FFT () is reflected around the Y-axis, so the negative half is the mirror of the positive half. This symmetry is caused by the input of real numbers (not complex numbers) into the transformation.

You can use this symmetry to make your Fourier Transform faster by computing only half of it. Scipy. FFT to RFFT ().

The nice thing about RFFT () is, it’s FFT (). Remember the previous FFT code:

yf = fft(normalized\_tone)  
xf = fftfreq(N, 1 / SAMPLE\_RATE)
Copy the code

Switch to RFFT () and the code remains basically the same, with a few key changes:

from scipy.fft import rfft, rfftfreq  
  
\# Note the extra 'r' at the front  
yf = rfft(normalized\_tone)  
xf = rfftfreq(N, 1 / SAMPLE\_RATE)  
  
plt.plot(xf, np.abs(yf))  
plt.show()
Copy the code

Since RFFT () only returns half the output FFT (), it uses a different function to get the frequency map, rfftfreq() rather than fftfreq().

RFFT () still produces complex output, so the code to draw its result remains the same. However, the graph should look like this because the negative frequency will disappear:

You can see that the figure above is just the front side of the spectrum generated by FFT (). RFFT () never computes the negative half of the spectrum, which makes it much better than using FFT ().

Usingrfft () can be twice as fast as using FFT (), but some input lengths are faster than others. If you know you can only use real numbers, then this is a speed trick worth knowing.

Now that you have the spectrum of the signal, you can proceed to filter it.

Filter the signal

One of the great things about the Fourier transform is that it’s reversible, so any changes you make to the signal in the frequency domain will apply when you transform it back into the time domain. You will use this to filter the audio and remove the high frequencies.

** Warning: ** The filtering techniques demonstrated in this section do not apply to real-world signals. For an explanation of why, see avoiding filter traps.

The returned value RFFT () represents the power of each frequency bin. If you set the power of a given bin to zero, the frequencies in that bin will no longer appear in the generated time domain signal.

Using the length Xf, the maximum frequency, and the fact that the frequency interval is evenly distributed, you can calculate the index of the target frequency:

\# The maximum frequency is half the sample rate  
points\_per\_freq = len(xf) / (SAMPLE\_RATE / 2)  
  
\# Our target frequency is 4000 Hz  
target\_idx = int(points\_per\_freq \* 4000)
Copy the code

You can then get rid of it by setting yf to 0 at the index around the target frequency:

yf\[target\_idx - 1 : target\_idx + 2\] = 0  
  
plt.plot(xf, np.abs(yf))  
plt.show()
Copy the code

Your code should generate the following diagram:

Since there’s only one peak, it seems to be working! Next, you will apply the Fourier Transform to return to the time domain.

Application of inverse FFT

Applying inverse FFT is similar to applying FFT:

from scipy.fft import irfft  
  
new\_sig = irfft(yf)  
  
plt.plot(new\_sig\[:1000\])  
plt.show()
Copy the code

Since you are using RFFT (), you need to use irfft() to apply the inverse. However, if you use FFT (), the inverse function will be ifft(). Your plot should now look something like this:

As you can see, you now have a sine wave oscillating at 400 Hz, and you have successfully eliminated 4000 Hz noise.

Again, you need to standardize the signal before writing it to a file. You can do what you did last time:

norm\_new\_sig = np.int16(new\_sig \* (32767 / new\_sig.max()))  
  
write("clean.wav", SAMPLE\_RATE, norm\_new\_sig)
Copy the code

When you listen to this file, you will hear the annoying noise disappear!

Avoid filter traps

The example above is more for educational purposes than practical use. Replicating the process on a real-world signal (a piece of music, for example) may introduce more buzz than it eliminates.

In the real world, you would use the filter design function in the Scipy. signal package to filter the signal. Filtering is a complex subject that involves a lot of math. See the Scientists and Engineers’ Guide to Digital Signal Processing for details.

Discrete cosine and sine transformations

Scipy.fft A tutorial on this module would be incomplete without an understanding of the discrete cosine Transform (DCT) and discrete Sine transform (DST). These two transformations are closely related to the Fourier transform, but operate entirely on real numbers. This means that they take one real-valued function as input and produce another real-valued function as output.

SciPy implements these transformations as DCT () and DST (). The I * and *n variants are inverse and n functional dimension versions, respectively.

DCT and DST are a bit like two halves of the Fourier Transform. This is not entirely true, as mathematics is much more complex, but it is a useful mental model.

So if DCT and DST are like half of Fourier Transform, why are they useful?

On the one hand, they are faster than a full Fourier Transform because they effectively do half the work. They can even be better than RFFT (). Best of all, they work entirely in real numbers, so you never have to worry about complex numbers.

Before you can learn how to choose between them, you need to understand even and odd functions. Even functions are symmetric about the y axis and odd functions are symmetric about the origin. To visualize this visually, check out this chart:

You can see that the even function is symmetric about the y axis. The odd function is symmetric with respect to y = -x, which is described as symmetric with respect to the origin.

When you compute the Fourier transform, you pretend that you are calculating its function to be infinite. The complete Fourier transform (DFT) assumes that the input function repeats indefinitely. However, DCT and DST assume that functions are extended symmetrically. DCT assumes that the function expands even symmetrically, while DST assumes that it expands odd symmetrically.

The following illustration illustrates how each transform imagines the function expanding to infinity:

In the figure above, the DFT repeats the function exactly as it is. DCT mirrors the function vertically to extend it, while DST mirrors it horizontally.

Note that the symmetry implicit in DST causes the function to jump around a lot. These are called ** discontinuities, ** and produce more high-frequency components in the resulting spectrum. Therefore, unless you know that your data has odd symmetry, you should use DCT instead of DST.

DCT is very common. There are more examples, but the JPEG, MP3, and WebM standards all use DCT.

conclusion

Fourier Transform is a powerful concept, used in a variety of fields, from pure mathematics to audio engineering and even finance. You are now familiar with the ** discrete Fourier Transform, ** and can use the scipy. FFT module nicely to apply it to filtering problems.

The last

Give the author a thumbs up if the article helps you

Continue to follow up with new Python articles, so don’t get lost