current position:Home>Use SciPy FFT for Fourier transform: Python signal processing

Use SciPy FFT for Fourier transform: Python signal processing

2022-02-02 06:51:41 Huawei cloud developer community

Abstract :Fourier transform It's a powerful concept , Used in various fields , From pure mathematics to audio engineering and even finance .

This article is shared from Huawei cloud community 《 Use scipy.fft Conduct Fourier Transform:Python signal processing 》, author : Yuchuan.

scipy.fft modular

Fourier transform is an important tool in many applications , Especially in scientific computing and Data Science . therefore ,SciPy Its implementation and related transformations have been provided for a long time . first ,SciPy Provided the scipy.fftpack modular , But then they updated their implementation and moved it to scipy.fft Module .

SciPy Full of functions . A more general introduction to the library , Please check out Scientific Python: Use SciPy To optimize .

install SciPy and Matplotlib

Before we start , You need to install SciPy and Matplotlib. You can do this in one of two ways :

  1. ** Use Anaconda install :** Download and install Anaconda Individual Edition. It has SciPy and Matplotlib, So once you follow the steps in the installer , You're done !

  2. ** Installation mode pip:** If you have pip install , Then you can install the library using the following command :

    $ python -m pip install -U scipy matplotlib

You can type... In the terminal python And run the following code to verify that the installation is valid :

>>>
>>> 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__.py
 Copy code 

This code imports SciPy and Matplotlib And print the location of the module . Your computer may display different paths , But as long as it prints the path , The installation is successful .

SciPy Now installed ! Now is the time to look scipy.fft The difference between and scipy.fftpack.

scipy.fft contrast scipy.fftpack

In the view SciPy When the document , You may encounter two modules that look very similar :

  1. scipy.fft

  2. scipy.fftpack

The scipy.fft The module is newer , Should take precedence over scipy.fftpack. You can go to SciPy 1.4.0 Read more about changes in the release notes for , But here's a quick summary :

  • scipy.fft There is an improved API.

  • scipy.fft Allow multiple worker, This can provide speed improvement in some cases .

  • scipy.fftpack Considered a legacy ,SciPy Suggest scipy.fft change to the use of sth. .

Unless you have good reason to use scipy.fftpack, Otherwise, you should stick to scipy.fft.

scipy.fft contrast numpy.fft

SciPy Fast Fourier transform (FFT) The implementation contains more functions , And ratio NumPy Your implementation is more likely to fix errors . If there is a choice , You should use SciPy Realization .

NumPy Maintained FFT Implement for backward compatibility , Although the author believes that functions such as Fourier transform are best placed in SciPy in . For more details , see also SciPy frequently asked questions .

Fourier Transform

Fourier Analysis is the study of how Mathematical functions Break down into a series of simpler The domain of trigonometric functions . Fourier transform is a tool in this field , Used to decompose a function into its component frequencies .

ok , This definition is very dense . For the purposes of this tutorial , Fourier transform is a tool , Allows you to get the signal and view the power of each frequency . Look at the important terms in this sentence :

  • One The signal It's information that changes over time . for example , Audio 、 Video and voltage traces are examples of signals .

  • nail frequency Is the rate at which something repeats . for example , The clock runs at one hertz (Hz) The frequency of ticking , Or repeat every second .

  • under these circumstances , power Only the intensity of each frequency .

The following figure is a visual demonstration of the frequency and power of some sine waves :

Said peak value High frequency of Sinusoids are closer together than those Low frequency Sine wave , Because they repeat more often . Described Low power The sine wave has a smaller peak than the other two sine waves .

To be more specific , Suppose you use Fourier transform for a recording of someone playing three notes on the piano at the same time . result spectrum Three peaks are displayed , One for each note . If the person plays one note softer than the others , Then the frequency intensity of this note will be lower than the other two .

This is the visual appearance of the piano example :

The highest note on the piano is quieter than the other two notes , Therefore, the spectrum of the note has a low peak .

Why Fourier Transform?

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

JPEG Compression uses a variant of Fourier transform to remove the high-frequency components of the image . Speech recognition uses Fourier transform and correlation transform to recover spoken language from the original audio .

Usually , If you need to check the frequency in the signal , Fourier transform is required . If it is difficult to process signals in the time domain , Then it is worth trying to use Fourier transform to move it to the frequency domain . In the next section , You will understand the difference between time domain and frequency domain .

Time Domain vs Frequency Domain

In the rest of this tutorial , You will see the term time domain and frequency domain. These two terms refer to two different ways of looking at signals , Or its component frequency , Or information that changes over time .

In the time domain , The signal is amplitude (y Axis ) Over time (x Axis ) Changing waves . You're probably used to seeing graphics in the time domain , For example, this :

Here are some audio images , It's a Time domain The signal . The horizontal axis represents time , The vertical axis represents the amplitude .

In frequency domain , The signal is represented as a series of frequencies (x Axis ), Each frequency has associated power (y Axis ). The following figure shows the above process Fourier Transform Audio signal after :

here , The previous audio signal is represented by its constituent frequency . Each frequency at the bottom has an associated power , Generate the spectrum you see .

More information about frequency domain , Please check out DeepAI Glossary entries .

Fourier Transform The type of

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

Even in this tutorial , You will often see DFT and FFT The two terms are used interchangeably . However , They are not exactly the same . Of ** fast fourier transform (FFT)** Is used to calculate the discrete Fourier transform (DFT) The algorithm of , and DFT It's the transformation itself .

You will be at scipy.fft Another difference seen in the library is the difference between different types of input .fft() Accept complex value input , and rfft() Accept real value input . Skip to using fast Fourier transform (FFT) Part to understand the complex and real numbers .

The other two transformations are related to DFT Is closely related to the : Discrete cosine transform (DCT) and Discrete sine transform (DST). You'll learn about these in the discrete cosine and Sine Transforms section .

Practical examples : Remove unwanted noise from audio

To help you understand the Fourier transform and what you can do with it , You will filter some audio . First , You will create an audio signal with a high hum , Then you'll use Fourier transform to remove the buzz .

Create signals

Sine waves are sometimes called ** pure tone ,** 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 use NumPy Directly generate . If you haven't used it before NumPy, So you can see what is NumPy?

This 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 code 

You will import and NumPy and Matplotlib, You can define two constants :

  1. SAMPLE_RATE Determines how many data points per second the signal uses to represent the sine wave . therefore , If the sampling rate of the signal is 10 Hz, And is 5 Second sine wave , Then it will have 10 * 5 = 50 The data points .

  2. DURATION Is the length of the generated sample .

Next , You define a function to generate a sine wave , Because you will use it many times in the future . This function uses frequency ,freq Then return to the... Used to draw the waveform x and y value .

Sine wave x The coordinates are 0 Evenly distributed between and DURATION, So the code uses NumPylinspace() To generate them . It needs a starting value 、 An end value and the number of samples to be generated . Set up endpoint=False It is important for the Fourier transform to work properly , Because it assumes that the signal is periodic .

np.sin() Calculate each x The value of the sine function at the coordinates . Multiply the result by the frequency , Make the sine wave oscillate at this frequency , Product times 2π To convert the input value to radians .

** Be careful :** If you haven't learned much trigonometry before , Or you need to review , So please check the trigonometry course of Khan College .

After defining the function , You can use it to generate a persistent 5 Of a second 2 Hertz sine wave , And use Matplotlib Draw it . Your sine wave diagram should look like this :

x The axis represents time in seconds , Because there are two peaks per second , You can see that the sine wave oscillates 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 , You will learn how to mix them .

Mixed audio signal

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

  1. Add the signals

  2. Standardization result

Before mixing the 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 code 

There is nothing new in this code example . It generates variables that are assigned to each other nice_tone and Midrange and treble noise_tone. You will use the treble as the noise you don't need , So it will multiply by 0.3 Reduce its power . Then the code adds these tones together . Please note that , You use underscores ( _) To discard by x The value returned generate_sine_wave().

The next step is to normalization , Or scale the signal to fit the target format . Because of how you will store audio later , Your target format is a 16 An integer , Range from -32768 To 32767:

normalized_tone = np.int16((mixed_tone / mixed_tone.max()) * 32767)

plt.plot(normalized_tone[:1000])
plt.show()
 Copy code 

ad locum , Code to zoom mixed_tone To make it suitable for 16 An integer , And then use NumPy Of np.int16. except mixed_tone Shrink it to... With its maximum value -1 Between and 1. When this signal is multiplied by when 32767, It's in -32767 Zoom between and 32767, Is roughly The scope of the np.int16. This code draws only the first 1000 sample , So that you can see the structure of the signal more clearly .

Your plot should be like this :

The signal looks like a distorted sine wave . The sine wave you see is generated by you 400 Hz tone , The distortion is 4000 Hz tone . If you look closely , You will find that the distortion is in the shape of a sine wave .

To listen to the audio , You need to store it in a format that the audio player can read . The easiest way is to use SciPy Of wavfile.write Method stores it in WAV In file .16 Bit integers are WAV The standard data type of the file , So you need to standardize the signal to 16 An integer :

from scipy.io.wavfile import write

# Remember SAMPLE_RATE = 44100 Hz is our playback rate
write("mysinewave.wav", SAMPLE_RATE, normalized_tone)
 Copy code 

This code will be written to mysinewave.wav You run Python Files in the directory of the script . then , You can use any audio player or even Python Listen to this file . You will hear lower tones and higher tones . These are your mixed 400 Hz and 4000 Hz Sine wave .

After completing this step , Your audio sample is ready . The next step is to use Fourier transform to remove the treble !

Use fast Fourier Transform (FFT)

It's time to use... On the generated audio FFT 了 .FFT It is an algorithm to realize Fourier transform , The spectrum of the signal in the time domain can be calculated , For example, 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 code 

This code will calculate the of the generated audio Fourier Transform And draw it . Before decomposing it , First look at the plot it produces :

You can see two peaks in the positive frequency and a mirror image of these peaks in the negative frequency . The positive frequency peak is located at 400 Hz and 4000 Hz, This corresponds to the frequency of the audio you input .

Fourier transform Your complex 、 A weak signal is converted into the frequency it contains . Because you only entered two frequencies , So only two frequencies came out . Negative positive symmetry is that Real value input Put in Fourier transform Side effects , But you'll hear more about it later .

In the first few lines , You can import scipy.fft The function to be used later , And define a variable N, Used to store the total number of samples in the signal .

After that is the most important part , Calculation Fourier transform

yf = fft(normalized_tone)
xf = fftfreq(N, 1 / SAMPLE_RATE)
 Copy code 

The code calls two very important functions :

  1. fft() Calculate the transformation itself .

  2. fftfreq() Calculation Each of the output of bin The frequency of the center fft(). Without this , You can't plot on the spectrum x Axis .

nail box Yes, it has been grouped , Like a histogram in a range of values . of bin Of For more information , See this signal processing stack swap problem . For the purposes of this tutorial , You can think of them as a single value .

Once you get Fourier transform The resulting value and its corresponding frequency , You can draw them :

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

The interesting part of this code is that you yf What you do before drawing . You call for np.abs() yes yf Because its value is complex .

nail The plural Is a number , It has two parts , namely real Department and virtual Ministry . There are many reasons why defining such numbers is useful , But what you need to know now is that they exist .

Mathematicians usually use a + bi Write the plural in the form of , among a It's the real part ,b It's the empty part .b hinder i Express b It's an imaginary number .

** Be careful :** Sometimes you will see using i Written plural , Sometimes you will see using j Written plural , for example 2 + 3 i and 2 + 3 j. The two are the same , but i Used more by mathematicians , and j Used more by engineers .

More about plural numbers , Please check the course of Khan college or the interesting page of Mathematics .

Because the plural has two parts , Plotting them against the frequency on the two-dimensional axis requires you to calculate a value from them . This is it. np.abs() Where I came in . It calculates the number of complex numbers √(a² + b²), This is the overall size of the two numbers , What matters is the individual values .

** Be careful :** By the way fft() A word of , You may have noticed , To be exact , The maximum frequency returned just exceeds 20,000 Hertz , namely 22050Hz. This value is exactly half of our sampling rate , It's called the Nyquist frequency .

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

Make it faster rfft()

fft() The output spectrum is wrapped around y Axial reflection , So the negative half is the mirror of the positive half . This symmetry is caused by the transformation input The set of real Numbers ( Not plural ) Caused by the .

You can use this symmetry , Make your... By counting only half of it Fourier transform faster .scipy.fft With rfft().

That's great rfft() yes , It is fft(). Remember the previous FFT Code :

yf = fft(normalized_tone)
xf = fftfreq(N, 1 / SAMPLE_RATE)
 Copy code 

Exchange in rfft(), The code remains basically the same , Just 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 code 

because rfft() Only half of the output is returned fft(), It uses different functions to get the frequency mapping ,rfftfreq() instead of fftfreq().

rfft() Still produce complex output , So the code that draws the result remains the same . however , The diagram shall be as follows , Because the negative frequency will disappear :

You can see that the above figure is just fft() The positive side of the generated spectrum .rfft() Never calculate the negative half of the spectrum , This makes it better than using fft().

usingrfft() Comparable using Twice as fast fft(), But some input lengths are faster than others . If you know you can only use real numbers , So this is a speed skill worth understanding .

Now you have the spectrum of the signal , You can continue to filter it .

Filter the signal

One of the great advantages of Fourier transform is that it is reversible , Therefore, any changes you make to the signal in the frequency domain will be applied when you transform it back to the time domain . You'll use this to filter audio and remove high frequencies .

** Warning :** The filtering techniques demonstrated in this section are not applicable to real-world signals . Explanation of the cause , See the section avoiding filter traps .

The value returned rfft() Represents the power of each frequency bin . If you will give bin The power of is set to zero , Then bin The frequency in will no longer appear in the generated time domain signal .

Use The length of xf、 The fact that the maximum frequency and frequency interval are 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 code 

then , You can set... At the index around the target frequency yf by 0 To get rid of it :

yf[target_idx - 1 : target_idx + 2] = 0

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

Your code should generate the following diagram :

Since there is only one peak , It seems to have worked ! Next , You will apply Fourier Transform Returns the time domain .

Apply inverse FFT

Apply inverse FFT Similar to application FFT:

from scipy.fft import irfft

new_sig = irfft(yf)

plt.plot(new_sig[:1000])
plt.show()
 Copy code 

Because you are using rfft(), You need to use irfft() To apply inverse . however , If you use fft(), Then the inverse function will be ifft(). Your plot should be like this now :

As you can see , You now have one with 400 Hz Oscillating sine wave , And you have successfully eliminated 4000 Hz noise .

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

norm_new_sig = np.int16(new_sig * (32767 / new_sig.max()))

write("clean.wav", SAMPLE_RATE, norm_new_sig)
 Copy code 

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

Avoid filtering traps

The above example is more for educational purposes , Instead of actually using . Signals in the real world ( For example, a piece of music ) Copying the process on may introduce more buzz than eliminating .

In the real world , You should use scipy.signal The filter in the packet is designed to filter the signal . Filtering is a complex subject involving a lot of Mathematics . For details , Check out the digital signal processing guide for scientists and engineers .

Discrete cosine and Sine Transforms

scipy.fft If you don't know the discrete cosine transform (DCT) And discrete sine transform (DST), The tutorial on this module will be incomplete . These two transformations are related to Fourier transform Is closely related to the , But it operates entirely on real numbers . This means that they take a real valued function as input , And generate another real valued function as the output .

SciPy Implement these transformations as dct() and dst(). Of i* and *n The variants are inverse and Ñ Functional dimension version of , , respectively, .

DCT and DST It's kind of like making up Fourier transform Two halves of . It's not exactly right , Because mathematics is much more complicated , But it's a useful mental model .

therefore , If DCT and DST It's like Fourier transform Half of , So why are they useful ?

One side , They are better than complete Fourier transform faster , Because they do half the work effectively . They can even be better than rfft(). most important of all , They work entirely in real numbers , So you never have to worry about the plural .

Before learning how to choose between them , You need to know even and odd functions . Even functions About y axial symmetry , and Odd functions About the origin symmetry . Imagine this intuitively , Please see the chart below :

You can see that even functions are about y axial symmetry . Odd functions about y = -x symmetry , This is described as About the origin symmetry .

When you calculate the Fourier transform , The function you pretend to be calculating is infinite . Complete Fourier transform (DFT) Suppose the input function repeats indefinitely . However ,DCT and DST Suppose the function is extended symmetrically .DCT Suppose that the function extends in even symmetry , and DST Suppose it extends with odd symmetry .

The following figure illustrates how each transformation imagines a function extending to infinity :

In the diagram above ,DFT The function is repeated as is .DCT Mirror the function vertically to extend it , and DST Mirror it horizontally .

Please note that ,DST Implicit symmetry can lead to large jumps in functions . These are called ** discontinuity ,** And generate more high-frequency components in the resulting spectrum . therefore , Unless you know that your data has odd symmetry , Otherwise you should use DCT instead of DST.

DCT Very often . There are more examples , but JPEG、MP3 and WebM All standards use DCT.

Conclusion

Fourier transform It's a powerful concept , Used in various fields , From pure mathematics to audio engineering and even finance . You are now familiar with ** discrete Fourier transform ,** And it can be used well scipy.fft The module applies it to filtering problems .

Click to follow , The first time to learn about Huawei's new cloud technology ~

copyright notice
author[Huawei cloud developer community],Please bring the original link to reprint, thank you.
https://en.pythonmana.com/2022/02/202202020651399059.html

Random recommended