Spectrum analysis of ADC sampling sequence based on stm32

   This article mainly introduces the right ADC The collected digital sequence FFT Spectrum analysis .

Determine sampling rate

   In addition to complying with Nyquist's sampling law, some problems need to be considered in determining the sampling rate . In a digital system , We can only do some finite discrete operations , For finite sequences , We can't do it DTFT, Can only do DFT. This requires Treat a finite sequence as a periodic sequence .

Normalized angular frequency

   The known sampling rate is $f_s$, Then a frequency is $f_0$（$f_0<f_s/2$） After the ideal cosine signal of is sampled, the sequence should be ：

$$x[n]=A\cos \left(2\pi f_0\cdot n{T}_0\right)=A\cos \left(2\pi f_0\cdot \frac{n}{f_s}\right)n\in \mathrm{Z}$$

   As can be seen from the formula above , The sequence we actually get is independent of the specific sampling rate or signal frequency , But by their ratio $f_0/f_s$ Decisive , It's like signal frequency $f_0$ Sampling frequency $f_s$ Normalized .
   When the signal is continuous , The frequency of the signal $f$ It can be any real number . And in this discrete case ,$f_0/f_s$ The frequency has been limited to 0~1 Between , The corresponding angular frequency is 0~$2\pi$.
   Combined with the periodicity of the spectrum, we can know $\pi$~$2\pi$ Spectrum and $-\pi$~$0$ The spectrum of is the same .

   For a real sequence ,0~$\pi$ The spectrum of corresponds to the frequency 0~$f_s/2$ Complex exponential sequence of , And the rest of the $\pi$~$2\pi$ The spectrum of is the complex conjugate of the former . The real sequence is formed by the superposition of two conjugate complex exponential sequences . According to this characteristic , In computing a real sequence DFT when , It's only half the length , Because the other half just takes complex conjugation .

Finite length sequence

   The previous analysis is all in “ Ideal ” Under the condition of ,“ Ideal ” It means that the cosine signal must always exist , From the beginning of the birth of the universe to the distant future, it always exists . Such a signal is not available in real life , The actual signal always has a beginning and an end .

Rectangular window

   Directly intercept several periodic points from the ideal sequence to form a finite length sequence , So the ideal sequence is multiplied by a rectangular window in the time domain , That is, their convolution in the frequency domain . The frequency domain representation of the ideal cosine sequence has been analyzed above , All that remains is to focus on the spectrum of the window function .

Calculate the rectangular window DTFT, Direct use DTFT Analytical formula of ：

$$X\left(e^{j\omega }\right)={\sum }_{n=0}^{N-1}{e}^{j\omega n}=\frac{1-e^{j\omega N}}{1-e^{j\omega }}$$

Then we can use the following formula to further simplify ,

$$1-e^{j2\theta }=1-\left(\cos \left(2\theta \right)+j\sin \left(2\theta \right)\right)=-j2\sin \left(\theta \right)e^{j\theta }$$

Finally get ：

$$X\left(e^{j\omega }\right)=e^{-j\omega (N-1)/2}\frac{\sin \left(\omega N/2\right)}{\sin \left(\omega /2\right)}$$

   stay $\omega =0$ The location of ,$X\left(e^{j\omega }\right)=1$, And in the $\omega$ yes $2\pi /N$ The integral times of ,$X\left(e^{j\omega }\right)=0$. So when the period of the signal is exactly $N$ When , There will be no spectrum leakage （ Please supplement the convolution process by yourself ）.

Rising cosine window （Hann window ）

   Rising cosine window is also called Hann window , It can be expressed as ：

$$w[n]=\frac{1}{2}\left(1-\cos \left(2\pi \frac{n}{N}\right)\right),n\in \left[0,N-1\right]$$

Calculate its DTFT：

$$X\left(e^{j\omega }\right)=\frac{1}{2}{\sum }_{n=0}^{N-1}\left(1-\cos \left(2\pi n/N\right)\right)e^{-j\omega n}$$

   I can't find the solution explicitly , But you can substitute some special values to observe . When $\omega =0$ when ,$X\left(e^{j\omega }\right)=1/2$, And when $\omega =\pm 2\pi /N$ when ,$X\left(e^{j\omega }\right)=1/4$.
  $X\left(e^{j\omega }\right)$ The result of convolution with the spectrum of ideal sequence will show spectrum leakage . For single frequency signals , Two peaks will be generated on both sides of the original single peak $1/2$ Peak amplitude peak , At the same time, the peak amplitude also attenuates to the original $1/2$.
   Although it is called “ Spectrum leakage ”, But it doesn't really leak out , But by convolution . Therefore, in general, adding up the leaked spectrum directly is only an approximation .
   Different window functions have different spectral characteristics , The side lobe of the rectangular window is relatively high , Once the spectrum leakage occurs, it will introduce a large error to the amplitude measurement of a specific frequency , But its frequency resolution is also the highest .Hann Although there will be spectrum leakage in the window , But its side lobe decays fast .

Sampling rate and sequence length

  DFT The minimum frequency that can be resolved by the obtained spectrum is called frequency resolution ,N Point sequence to do DFT Can get N Spectrum of points , Therefore, the frequency resolution is $f_s/N$. It is best to choose the appropriate sampling rate and sequence length so that the frequency of the signal to be measured is exactly an integral multiple of the frequency resolution . If the frequency of the signal to be measured is exactly an integral multiple of the frequency resolution , Then you can use rectangular windows directly ; conversely , If the frequency of the signal to be measured is not an integral multiple of frequency resolution , Then you need to add windows according to your needs .

Analog front end

   In determining ADC After the sampling rate of ,ADC The cut-off frequency of the previous anti aliasing filter should also be set at $f_s/2$ within . But it should also be greater than the frequency of the signal you care about .
   This oversampling method is generally used to sample some low-frequency signals ,ADC The previous anti aliasing filter is a low-pass filter , The sampling rate is greater than twice the maximum frequency of the signal . But the sampling theorem tells us that as long as the sampling rate is greater than Signal bandwidth The sampled signal can be recovered by twice . So there is also a bandpass system ,ADC The former is a bandpass filter , Using the characteristics of spectrum periodic extension of digital sequence , Use low-frequency sampling rate to collect high-frequency signals .

Simulated signal

   stay ADC It hasn't been adjusted yet , Or sometimes I want to generate a single frequency simulation sequence by myself , You can write directly at the sampling rate $f_s$ Next , frequency $f_0$ Cosine sequence of ：

$$x[n]=\cos \left(2\pi \frac{f_0}{f_s}n\right)$$

DSP Library function

   In the previous one article in , I mentioned using CubeMX After project generation , stay Drivers/CMSIS/Lib It's in the catalog DSP The library files , So I don't need to import the library from the outside DSP The library .
   But these two library files are a little different , I am using CFFT I found that only the library files imported from the outside can be calculated correctly . If you use the existing library file , The program is in link Then it seems that something strange has covered the entry of the reset address , Come up here hardfault. Here, so it's better to check this again DSP library .

  CMSIS Of DSP The library provides two kinds of FFT Function of , The set of real Numbers FFT And the plural FFT, Both have their own characteristics . Can be used to calculate by ADC Of the sampled sequence FFT.

The set of real Numbers FFT

   The set of real Numbers FFT Fast calculation , It's our first choice ！ It can handle a length of 32, 64, 128, 256, 512, 1024, 2048, 4096 The real number FFT.
   The input sequence is a continuous real number , I thought the input was also to be plural FFT equally , The real part and the imaginary part are stored at intervals . But in fact, what it needs is continuous real numbers .
   Not the same address operation , Input and output need to be stored in different places . The output sequence is a complex sequence , The length is half of the input sequence , Because the results of the other half and the first half are conjugate , There is no need to double count .
   Because the complex number is stored in real part and imaginary part at intervals , A complex number takes up twice the storage space of a real number . So the actual physical space occupation of input and output is the same .
   For example 1024 Dot floating point RFFT, The input is a piece of length 1024 The first address of the address space of floating-point numbers , The output also needs a space of the same size , But the storage will be 512 Plural .
   Output storage is not all plural , A sequence of real numbers $x[n]$ do FFT And then get $X[k]$,$X[0]$ It must be a real number ,$X[N/2]$ It must also be a real number , These two real numbers exist together in the position of the first complex number .
   Sample code ：

arm_rfft_fast_instance_f32 xInstFFT;
arm_rfft_fast_init_f32(&xInstFFT, 1024);
arm_rfft_fast_f32(&xInstFFT, pSrc, pDst, 0);

   Calling RFFT Before , You need to initialize an instance first , It mainly sets the lookup table of rotation factor .RFFT All the time bit reverse, In this way, the input and output are in the order we are used to , and CFFT Yes, you can choose whether to bit reverse.arm_rfft_fast_f32 The last parameter in is choose to FFT still IFFT.

The plural FFT

#include "arm_const_structs.h"
arm_cfft_f32(arm_cfft_sR_f32_len1024, pData, 0, 1);

   The plural FFT The data required to be input is stored alternately in real part and imaginary part , And it is the same address operation , The output results will overwrite the input data . Compared with RFFT, It doesn't need to be initialized , It only needs include arm_const_structs.h This header file , You can find ready-made examples .

   This header file is in CMSIS Of DSP/Include Under the path , Although we can't use what is provided here DSP The library files , But the header file can still be used . The penultimate parameter is selection FFT still IFFT; The last parameter is whether bit reverse.

Hann window

  Hann The coefficient of the window can be used Matlab Calculation , These coefficients can be stored as constants , It is convenient to call directly when adding windows .

   open Matlab Window designer , choice Hann window , Then select “ periodic ”. Save the coefficients to the workspace , Then find a way to write C/C++ Just go to the source file of .

function [w] = genHann()
    w = hann(1024, 'periodic');
end

   I use Matlab Coder Turn the above function into C Code , But it's also troublesome , It's better to copy it directly from the workspace and write it to the source file .