Direct inverse discrete Fourier transform. Discrete Fourier transform

This is one of the Fourier transforms widely used in digital signal processing algorithms (its modifications are used in audio compression in MP3, image compression in JPEG, etc.), as well as in other areas related to the analysis of frequencies in discrete (for example, digitized analog ) signal. Discrete conversion Fourier requires as input discrete function. Such functions are often created by sampling (sampling values ​​from continuous functions). Discrete Fourier transforms help solve partial differential equations and perform operations such as convolutions. Discrete Fourier transforms are also actively used in statistics, in the analysis of time series. Transformations can be one-dimensional, two-dimensional and even three-dimensional.

Direct conversion:

Reverse conversion:

Designations:

§ N- the number of signal values ​​measured over a period, as well as the number of decomposition components;

§ - measured signal values ​​(at discrete time points with numbers , which are input data for direct conversion and output data for reverse conversion;

§ - N complex amplitudes of sinusoidal signals composing the original signal; are output data for direct conversion and input data for reverse conversion; since amplitudes are complex, it is possible to calculate both amplitude and phase from them;

§ is the usual (real) amplitude of the kth sinusoidal signal;

§ arg( X k) - phase of the kth sinusoidal signal (argument of a complex number);

§ k- frequency of the kth signal, equal to , where T- the period of time during which the input data was taken.

From the latter it is clear that the transformation decomposes the signal into sinusoidal components (which are called harmonics) with frequencies from N oscillations per period to one oscillation per period. Since the sampling frequency itself is equal to N samples per period, high-frequency components cannot be displayed correctly - a moiré effect occurs. This leads to the fact that the second half of the N complex amplitudes is, in fact, a mirror image of the first and does not carry additional information.

Let's consider some periodic signal x(t) with a period equal to T. Let us expand it into a Fourier series:

Let us sample the signal so that there are N samples per period. Let's represent the discrete signal in the form of samples: x n = x(tn), where , then these readings through the Fourier series will be written as follows:

Using the relation: , we get:

Where

So we got inverse discrete Fourier transform.

Let us now scalarly multiply the expression for x n on and we get:


Here we use: a) expression for the sum finite number members (exhibitor) geometric progression, and b) the expression of the Kronecker symbol as the limit of the ratio of the Euler functions for complex numbers. It follows that:

This formula describes direct discrete Fourier transform.

In the literature, it is customary to write the multiplier in the inverse transformation, and therefore the transformation formulas are usually written in the following form:

The discrete Fourier transform is a linear transformation that transforms a vector of time samples into a vector of spectral samples of the same length. So the transformation can be implemented as a multiplication square matrix to vector:

The program code for direct and inverse Fourier transform is given. The fast Fourier transform is considered.

Discrete Fourier Transform (DFT) is powerful tool analysis, which is widely used in the field of digital signal processing (DSP). There are direct and inverse Fourier transforms. The direct discrete Fourier transform converts a signal from the time domain to the frequency domain and is used to analyze the frequency spectrum of the signal. The inverse transformation does exactly the opposite: by frequency spectrum signal reconstructs the signal in the time domain.

To calculate the Fourier transform, an accelerated calculation procedure is usually used - the so-called. fast Fourier transform (FFT). This allows you to significantly reduce processor time for fairly complex and resource-intensive mathematical calculations.

1 Complex numbers

First, we need a helper class that will describe complex numbers. Complex numbers are special kind numbers in mathematics. Every complex number consists of two parts - real and imaginary. Now it is enough for us to know about complex numbers in relation to DFT that the real part of a complex number stores information about the signal amplitude, and the imaginary part stores information about the phase.

Class code for describing complex numbers(turns around) """ """ Complex number. """ Public Class ComplexNumber """ """ Real part of a complex number. """ Public Real As Double = 0 """ """ Imaginary part of a complex number. """ Public Imaginary As Double = 0 Public Sub New() Real = 0 Imaginary = 0 End Sub """ """ Creates a complex number. """ """ Real part of a complex number. """ Imaginary part of a complex number. Public Sub New(ByVal r As Double, Optional ByVal im As Double = 0) Real = r Imaginary = im End Sub Private usCult As New Globalization.CultureInfo("en-US") "we use the culture "en-US" so that the whole and Fractional parts were separated by a dot, not a comma """ """ Returns a string consisting of a real and an imaginary part, separated by a tab character. """ Public Overrides Function ToString() As String Return (Real.ToString(usCult) & ControlChars.Tab & Imaginary.ToString(usCult)) End Function End Class

2 Direct discrete fast Fourier transform

An array of complex numbers is passed to the function input. The real part of which represents an arbitrary discrete signal, with readings at regular intervals. The imaginary part contains zeros. The number of samples in the signal must be equal to a power of two. If your signal is shorter, then pad it with zeros to a multiple of 2: 256, 512, 1024, etc. The longer the signal, the higher the frequency resolution of the calculated spectrum.

Code for calculating direct fast Fourier transform in VB.NET(turns around) """ """ Calculates the spectrum of a signal using the Fast Fourier Transform method. Use only (N/2+1) return values ​​(up to half the sample rate). """ """ A signal containing a number of samples that is a multiple of a power of two and consists of a real and an imaginary part. All imaginary parts of the signal are filled with zeros. """ Returns an array of complex spectrum numbers. """ Only the first N/2+1 are significant, the rest are the symmetrical part corresponding to negative frequencies. """ The first value of the spectrum is the constant component, the last one corresponds to half the sampling frequency (Nyquist frequency). """ Values ​​above half the sampling frequency - do not use. """ Public Shared Function FFT(ByVal signal As ComplexNumber()) As ComplexNumber() Dim order As Integer = signal.Length "DFT order CheckFftOrder(order) "Check that the order equal to the power twos Dim spectrumLen As Integer = order \ 2 Dim j As Integer = spectrumLen "Bit-reverse sort: For i As Integer = 1 To order - 2 If (i< j) Then Dim tmpRe As Double = signal(j).Real Dim tmpIm As Double = signal(j).Imaginary signal(j).Real = signal(i).Real signal(j).Imaginary = signal(i).Imaginary signal(i).Real = tmpRe signal(i).Imaginary = tmpIm End If Dim k As Integer = spectrumLen Do Until (k >j) j -= k k \= 2 Loop j += k Next "Loop through expansion levels: For level As Integer = 1 To CInt(Math.Log(order) / Math.Log(2)) Dim lvl As Integer = CInt (2 ^ level) Dim lvl2 As Integer = lvl \ 2 Dim tmp As Double = Math.PI / lvl2 Dim sr As Double = Math.Cos(tmp) Dim si As Double = -Math.Sin(tmp) Dim tr As Double = 0 Dim ur As Double = 1 Dim ui As Double = 0 For jj As Integer = 1 To lvl2 "Cycle through spectra within a level For i As Integer = (jj - 1) To (order - 1) Step lvl "Cycle through individual "butterflies" Dim ip As Integer = i + lvl2 tr = signal(ip).Real * ur - signal(ip).Imaginary * ui "Butterfly operation" Dim ti As Double = signal(ip).Real * ui + signal (ip).Imaginary * ur signal(ip).Real = signal(i).Real - tr signal(ip).Imaginary = signal(i).Imaginary - ti signal(i).Real = signal(i).Real + tr signal(i).Imaginary = signal(i).Imaginary + ti Next tr = ur ur = tr * sr - ui * si ui = tr * si + ui * sr Next Next "Fill in the array of complex numbers processed by the FFT: Dim spectrum(order - 1) As ComplexNumber For i As Integer = 0 To order - 1 With signal(i) spectrum(i) = New ComplexNumber(.Real, .Imaginary) End With Next Return spectrum End Function

3 Inverse discrete fast Fourier transform

The inverse discrete Fourier transform (IDFT), one of the calculation stages, includes a direct DFT on an array of complex numbers, where the imaginary part is the inversion with respect to the X-axis of the imaginary part of the spectrum.

Code for calculating the inverse fast Fourier transform in VB.NET(turns around) """ """ Restores a signal from its spectrum using the inverse fast Fourier transform method. """ """ A signal spectrum containing a number of samples that is a multiple of a power of two and consists of a real and an imaginary part. Public Shared Function InverseFFT(ByVal spectrum As ComplexNumber()) As ComplexNumber() Dim order As Integer = spectrum.Length "Inverse DFT order. CheckFftOrder(order) "Changing the arithmetic sign of the elements of the imaginary part: For i As Integer = 0 To spectrum. Length - 1 spectrum(i).Imaginary = -spectrum(i).Imaginary Next "Direct FFT calculation: Dim directFFT As ComplexNumber() = FFT(spectrum) "Division by order in the time domain with changing the arithmetic sign of the imaginary part: Dim signal (directFFT.Length - 1) As ComplexNumber For i As Integer = 0 To directFFT.Length - 1 Dim ReX As Double = directFFT(i).Real / order Dim ImX As Double = -directFFT(i).Imaginary / order signal( i) = New ComplexNumber(ReX, ImX) Next Return signal End Function

And of course, we will describe the method used, which checks the number of elements of the passed array:

"""

""" Checks whether the FFT order is a power of two, and if not, throws an exception. """ """ FFT order. Private Shared Sub CheckFftOrder(ByVal order As Integer) Dim chk As Double = Math.Abs(Math.Floor(Math.Log(order, 2)) - Math.Log(order, 2)) If (chk > 0.0001) Then Throw New ArgumentException(String.Format("The length of the array ((0)) is not a power of two.", order)) End If End Sub

4 Checking forward and reverse Fourier transform

Now let's check that our functions work. To do this, let's pass an arbitrary signal through the direct Fourier transform mechanism, and then “assemble” it back using the inverse Fourier transform. The reconstructed signal should practically coincide with the original one. Rounding errors that occur when working with numbers on a computer do occur, so the signals will not be completely identical, but their deviation from each other should be negligible.

For example, let’s take the sine function as the source signal and generate data with a length of 128 samples like this:

Dim cn(127) As ComplexNumber For i As Integer = 0 To cn.Length - 1 cn(i) = New ComplexNumber(Math.Sin(i * 3 * Math.PI / 180)) Next

We get this signal:

Here, the X axis is the number of samples in the time domain, and the Y axis is the amplitude. Please note that the signal consists only of real parts, and the imaginary part throughout the entire segment is equal to “0”.

Now let's pass this signal to the input of the FFT() function. Using the arrays of complex numbers obtained during the direct Fourier transform, we will construct two graphs - the real (Re) and imaginary (Im) parts of the spectrum:


Here along the X axis are readings in frequency domain, along the Y axis - amplitude. To obtain the actual frequency values, it is necessary to calculate them, taking into account that the "0" of the Y axis corresponds to zero frequency, the maximum of the Y axis corresponds to the sampling frequency.

We will transfer the resulting signal spectrum to the inverse Fourier transform function IFFT(). Let's get an array of complex numbers, where the real part will contain the reconstructed signal:


As you can see, the reconstructed signal completely repeats the original one.

It is convenient to analyze many signals by decomposing them into sinusoids (harmonics). There are several reasons for this. For example, the human ear works in a similar way. It decomposes sound into individual vibrations of different frequencies. Additionally, it can be shown that sinusoids are " own functions"linear systems (since they pass through linear systems, without changing the shape, but can only change the phase and amplitude). Another reason is that Kotelnikov's theorem is formulated in terms of the signal spectrum.

The Fourier transform is the decomposition of functions into sinusoids (hereinafter we also call cosine functions sinusoids, since they differ from “real” sinusoids only in phase). There are several types of Fourier transform.

1. Non-periodic continuous signal can be expanded into a Fourier integral.

2. A periodic continuous signal can be expanded into an infinite Fourier series.

3. A non-periodic discrete signal can be expanded into a Fourier integral.

4. A periodic discrete signal can be expanded into a finite Fourier series.

A computer can only work with a limited amount of data, therefore, it can really only calculate the last type of Fourier transform. Let's take a closer look at it.

Complex DFT

So far we have considered DFTs from real signals. Let us now generalize the DFT to the case of complex signals. Let x[n], n=0,…,N-1 - the original complex signal consisting of N complex numbers. Let us denote X[k], k=0,…N-1 - its complex spectrum, also consisting of N complex numbers. Then the following formulas for direct and inverse transformations Fourier:

If we decompose a real signal into a spectrum using these formulas, then the first N/2+1 complex coefficients of the spectrum will coincide with the spectrum of the “usual” real DFT, presented in “complex” form, and the remaining coefficients will be their symmetric reflection relative to half the sampling frequency. For cosine coefficients the reflection is even, and for sine coefficients it is odd.

2D DFT

For images that are a two-dimensional signal, the spectrum is also a two-dimensional signal. The basis functions of the Fourier transform have the form:

Moreover, the phases can also be different. In the image, each of these basis functions represents a wave of a certain frequency, a certain orientation, and a certain phase.

Here N 1 xN 2 is the size of the original signal, which is also the size of the spectrum. k 1 and k 2 are the numbers of the basis functions (the numbers of the coefficients of the two-dimensional DFT at which these functions are found). Since the spectrum size equal to size source signal, then k 1 = 0,…,N 1 -1; k 2 = 0,…,N 2 -1.

n 1 and n 2 are variable arguments of the basis functions. Since the domain of definition of the basis functions coincides with the domain of definition of the signal, then n 1 = 0,...,N 1 -1; n 2 = 0,…,N 2 -1.

2D DFT (in complex form) is determined the following formulas(here x is the original signal and X is its spectrum):

Direct calculation of a two-dimensional DFT using the above formulas requires enormous computational costs. However, it can be proven that the two-dimensional DFT has the separability property, i.e. it can be calculated sequentially from two dimensions.

To calculate a two-dimensional DFT, it is enough to calculate the one-dimensional complex DFTs of all rows of the image, and then calculate the one-dimensional complex DFTs of all columns in the resulting “image”.

In this case, the results of all one-dimensional complex DFTs must be written in place of the original data for these DFTs. For example, when calculating the one-dimensional DFT of the first row of an image, you need to write the DFT result in the first row of this image (it has the same size). To do this, you need to store each “pixel” as a complex number.

Thus, efficient algorithm Calculating the DFT of an image consists of calculating one-dimensional FFTs first from all rows, and then from all columns of the image.

Let f(x 1 , x 2) – a function of two variables. By analogy with the one-dimensional Fourier transform, we can introduce a two-dimensional Fourier transform:

The function for fixed values ​​of ω 1, ω 2 describes plane wave in the plane x 1 , x 2 (Figure 19.1).

The quantities ω 1, ω 2 have the meaning of spatial frequencies and dimension mm−1, and the function F(ω 1, ω 2) determines the spectrum of spatial frequencies. A spherical lens is capable of calculating the spectrum of an optical signal (Figure 19.2). In Figure 19.2 the following notations are introduced: φ - focal length,

Figure 19.1 - To determine spatial frequencies

The two-dimensional Fourier transform has all the properties of the one-dimensional transform, in addition, we note two additional properties, the proof of which easily follows from the definition of the two-dimensional Fourier transform.


Figure 19.2 – Calculation of the spectrum of an optical signal using
spherical lens

Factorization. If a two-dimensional signal is factorized,

then its spectrum is also factorized:

Radial symmetry . If the two-dimensional signal is radially symmetric, that is

Where is the zero-order Bessel function. The formula that defines the relationship between a radially symmetric two-dimensional signal and its spatial spectrum is called the Hankel transform.


LECTURE 20. Discrete Fourier transform. Low Pass Filter

The direct two-dimensional discrete Fourier transform (DFT) transforms an image given a spatial coordinate system (x, y), into a two-dimensional discrete image transformation specified in a frequency coordinate system ( u,v):

The inverse discrete Fourier transform (IDFT) has the form:

It can be seen that the DFT is complex transformation. The modulus of this transformation represents the amplitude of the image spectrum and is calculated as the square root of the sum of the squares of the real and imaginary parts of the DFT. Phase (phase shift angle) is defined as the arctangent of the ratio of the imaginary part of the DFT to the real part. Energy spectrum equal to square amplitude of the spectrum, or the sum of squares of the imaginary and real parts of the spectrum.



Convolution theorem

According to the convolution theorem, the convolution of two functions in the spatial domain can be obtained by the ODFT of the product of their DFT, that is

Filtering in the frequency domain allows you to use the DFT of the image to select the frequency response of the filter that provides the necessary image transformation. Let's look at the frequency characteristics of the most common filters.

In radio engineering, the concept of convolution of two signals is often used. For example, the signal at the output of a four-port network can be found by convolving the input signal and the impulse response of the four-port network. Since discrete and digital signals were considered, we define the concept convolution for discrete signals , or discrete convolution.

Let there be a discrete signal x D (t), consisting of N counts x k, and a discrete signal y d (G), consisting of N counts u k, Then discrete convolution these two signals are called signal zA(t), for which

Discrete signals have become widespread in the creation of pulse modulation systems.

In the simplest case, the sampling device is a gated cascade (key) that opens for a while t and with period A (Fig. 4.7).


Rice. 4.

The sampling interval A can be constant (uniform sampling) or variable (adaptive sampling). The most common form of discretization is uniform, which is based on Kotelnikov’s theorem.

Pulse modulator - This is a device with two inputs, one of which is supplied analog signal, and the second receives short synchronizing pulses with a repetition period A. In this case, at the moment the synchronizing pulse arrives, the instantaneous value of the signal hp(g) is measured. A sequence of pulses appears at the output of the modulator, each of which has an area proportional to the corresponding reference value of the analog signal (Fig. 4.7).

Signal Hmpn ( t) at the output of the pulse modulator is called modulated pulse sequence(MIP). Mathematically, the MIP is written as follows:

A MIP spectral density expressed through spectral density analog signal as follows:

The discrete signal model assumes that sample values ​​of an analog signal can be obtained at an unlimited number of points on the time axis. In practice, processing is always carried out over a finite time interval.

Let us consider the features of the spectral representation of a discrete signal, defined on an interval by its own samples x 0 ,x x ,...,x N _ x . Full number counts N - T/ A.

The technique for studying such discrete signals is that the resulting sample of reference values ​​is mentally repeated infinite number once. As a result, the signal becomes periodic (Fig. 4.8).

By matching such a signal mathematical model, you can use the Fourier series expansion and find the corresponding amplitude coefficients. The combination of these coefficients forms the spectrum of a discrete periodic signal.


Rice. 4.8.

Let us write the model of a limited periodic signal in the form of a sequence of delta pulses:

Let us expand the signal Xmip (0 into a Fourier series:

Is this a change of variables? = f / A. Finally we get

This formula determines the sequence of coefficients forming discrete Fourier transform (DFT) of the signal under consideration.

DFT has the following properties:

1. DFT is a linear transformation, i.e. if z k = a x k + /? u k, That

C "Z P ~ ^ C X p R Su p .

2. Number of different coefficients Cq,Ci,...,C n _i is equal to the number N counts for the period, with n = N coefficient C N= C 0 .

3. C 0 is the average value of all readings C 0 = - ^x To.

N to=o

  • 4. If N- even number, That With N = -^(-1) k x k.
  • 7 ^ ?=o
  • 5. If the readings x k - real numbers And N is an even number, then C N = C* N, / = 0; L/7 2 -1.
  • -+i - -i
  • 6. If y k =x k+m , m = l;JV-l,TO C, t =C, * e ~ j2rrkm,N .
  • 2 tf-l
  • 7. If z k= -> T0 C z k =C X k C y k

iy/i=0

DFT is used to calculate the spectra of functions specified by tables or graphs, process experimental data, find the signal at the output of a discrete filter, etc.

If based on readings x 0 ,x l ,...,x N _ l DFT coefficients are found for some signal C 0 ,Ci,... 9 C n/2 , then using them you can reconstruct an analog signal with a limited spectrum x(t). The Fourier series of such a signal has the form (for even N)

where |Q| - module of DFT coefficients; =arg - phase angle (argument)

DFT coefficients. First harmonic frequency: f= -/ in = - = -/i- odd N the last term in formula (4.17) is equal to:

To calculate discrete samples x k Based on the available DFT coefficients, there is the following formula:

This formula is called inverse discrete Fourier transform (IDFT).



Did you like the article? Share with your friends!