current position：Home>Use SciPy FFT for Fourier transform: Python signal processing
Use SciPy FFT for Fourier transform: Python signal processing
20220202 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 ：

** 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 ！

** 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/distpackages/scipy/__init__.py
>>> print(matplotlib.__file__)
/usr/local/lib/python3.6/distpackages/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 ：

scipy.fft

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 highfrequency 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 ：

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 .

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 ：

Add the signals

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 ：

fft() Calculate the transformation itself .

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 twodimensional 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 realworld 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 highfrequency 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
The sidebar is recommended
 [Python] numpy notes
 24 useful Python tips
 Pandas table beauty skills
 Python tiktok character video, CV2 module, Python implementation
 I used Python to climb my wechat friends. They are like this
 20000 words take you into the python crawler requests library, the most complete in history!!
 Answer 2: why can you delete the table but not update the data with the same Python code
 [pandas learning notes 02]  advanced usage of data processing
 How to implement association rule algorithm? Python code and powerbi visualization are explained to you in detail (Part 2  actual combat)
 Python adds list element append() method, extend() method and insert() method [details]
guess what you like

python wsgi

Introduction to Python gunicorn

Python dictionary query key value pair methods and examples

Opencv Python reads video, processes and saves it frame by frame

Python learning process and bug

Imitate the up master and realize a live broadcast room controlled by barrage with Python!

Essence! Configuration skills of 12 pandas

[Python automated operation and maintenance road] path inventory

Daily automatic health punch in (Python + Tencent cloud server)

[Python] variables, comments, basic data types
Random recommended
 Spring boot calls Python interface
 Using Python to make a key recorder
 Python combat case, pyGame module, python implementation routine confession artifact vs no routine confession artifact
 Python series tutorial 132  why use indentation syntax
 10 minutes to learn how to play excel easily with Python
 Python develops a color dynamic twodimensional code generator in one hour, and uses the virtual environment to package and release the EXE program.
 Elimination of grammar left recursion in Python
 Python testing  the patches in Python
 Python image processing, CV2 module, OpenCV to achieve target tracking
 How to send alarm notification to nail in Python?
 Introduction to pandas operation
 Mail sending, SMTP and exchange sending in Python 3
 Show your hand and use Python to analyze house prices
 The strongest Python visualization artifact, none of them
 8 practical Python skills that are easy to use and don't have to suffer a loss for half a year
 Tips: teach you to generate 30 cool dynamic interactive charts with one click of pandas
 I use one line of Python code to dynamically load dependencies
 Blow up this pandas GUI artifact and automatically turn the code!
 Getting started exploring and analyzing data using Python
 Python image processing, CV2 module, OpenCV to achieve template matching
 Teach you three Python Programming Tips
 Python code reading (Chapter 69): case conversion of initial letters
 Python tutorial series 133  several special grammars
 Dry goods  Python operation: Python controls Excel to realize office automation
 Understand the principle of affine transformation and python implementation code
 A little trick every day, python can easily convert PDF to text and bid farewell to copy and paste
 Climb Conan barrage with Python + gephi to sort out the main plot
 As a programmer, do you know Python variable reference
 One line of code, take you to learn Python
 Summary of 22 advanced Python knowledge points, dry goods!
 [Python learning] nanny level teaching parsing and parsing XML in Python
 Use of Python JSON module
 "Python" guide to using itertools of Python standard library
 What are the functions and classifications of Python interpreters
 Python implements four schemes for timed tasks, lazy artifact
 Python can actually realize the freedom of punch in?
 What about Python memory leak? Pit filling troubleshooting tips
 Object oriented programming in Python
 Quick start  Python Basics
 [target detection (8)] a thorough understanding of the loss function of the regression box of target detection  the principles of IOU, giou, Diou and ciou and Python code