A high-precision fast implementation method, system and device for fractional fourier transform of digital signal processing

Abstract

A high-precision fast implementation method, system and device for fractional Fourier transform of digital signal processing, the method comprising: determining the analysis range of the angle of fractional Fourier transform, the total size of the time delay analysis range and the time delay analysis step according to the situation and requirements of the signal to be processed; calculating the extended analysis point number and the point number expansion multiple according to the angle range, the total size of the time delay analysis range and the time delay analysis step, and zero-padding the discretized input signal; calculating the digital angular frequency resolution interval and the index value; calculating the synthesized signal according to the digital angular frequency resolution interval, the index value and the input signal; performing anti-aliasing filtering and downsampling on the synthesized signal to obtain the downsampled signal; performing fast Fourier transform on the downsampled signal to obtain the preliminary result; performing inversion, downsampling and truncation on the preliminary result to obtain the final fractional Fourier transform result; the system and device are used to implement the method; the present application has the characteristics of high calculation precision, fast calculation, convenient and flexible use.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of radar and underwater acoustic signal processing, and particularly relates to a high-precision fast implementation method, system and equipment for fractional Fourier transform of digital signal processing. BACKGROUND

[0002] As a popular signal time-frequency domain analysis method, fractional Fourier transform (FrFT) is widely used in the field of signal processing, and it is particularly suitable for analyzing linear frequency modulation signals. By performing fractional Fourier transform on a signal, the signal is projected into an alpha-u two-dimensional transform domain, where alpha represents the projection angle, which directly matches the frequency modulation slope of the signal, and u is related to the time delay of the signal. Through fractional Fourier transform, the time-frequency distribution of the signal can be easily analyzed, which facilitates subsequent processing such as signal detection and recognition.

[0003] If the fractional Fourier transform is calculated directly according to the definition, it will have a large complexity, which is not conducive to practical engineering application. In order to facilitate the practical application of fractional Fourier transform, various fast calculation techniques and systems of this method have been derived. The main principle of these fast calculation techniques is as follows: for each angle alpha, the corresponding fractional Fourier transform is decomposed into two or more linear convolutions of signals, and then the linear convolution operation is realized by fast Fourier transform (FFT) technology to realize fast calculation of fractional Fourier transform. The main difference between the above different techniques lies in the different ways of decomposing the fractional Fourier transform, but the same point mainly reflects the decomposition into linear convolution form of different signal components.

[0004] In the FFT fast calculation process of linear convolution, multiple FFT and IFFT (inverse fast Fourier transform) operations are required. Specifically, each convolution term is generally transformed into the discrete Fourier frequency domain by FFT, the transformed terms are multiplied in the discrete Fourier frequency domain, and then the multiplication result is transformed back to the time domain by IFFT, thereby obtaining the fast calculation result of linear convolution. Although FFT and IFFT techniques have high calculation efficiency, multiple applications of FFT and IFFT will obviously increase the calculation amount by several times compared with single application. On the other hand, in the existing FrFT fast calculation techniques, the two or more signals (or functions) decomposed into linear convolution are mainly in the form of linear frequency modulation. In the field of wideband or ultra-wideband signal processing, due to the characteristics of wide frequency band and long signal sequence, the decomposition of signals or functions and their operations will inevitably bring about a non-negligible calculation error, which will have a significant impact on the result, resulting in loss of calculation precision.

[0005] The Chinese patent application with the publication number CN 107644004A discloses a digital signal processing method and device based on a fast calculation method of discrete fractional Fourier transform. However, since the method needs to calculate the eigenvalues and eigenvectors of the discrete Fourier transform kernel matrix, the calculation complexity is increased, and the calculation efficiency is lost.

[0006] The Chinese patent application with the publication number CN 105783974A discloses a method and system for detecting and estimating parameters of a linear frequency modulation signal. The method realizes the calculation of fractional Fourier transform through rotation transform and fast Fourier transform (FFT). However, the resolution accuracy and the transform domain range of the corresponding transform domain calculation result are strictly limited by the sampling point number N and the sampling frequency fs of the signal, the design flexibility is lost, and the method performs FFT on the entire sampling frequency range fs, which increases the calculation redundancy and also loses the calculation efficiency. SUMMARY

[0007] In order to overcome the above-mentioned shortcomings of the prior art, the purpose of the present application is to provide a fractional Fourier transform high-precision fast implementation method, system and device for digital signal processing. The method realizes faster calculation compared with the prior art by newly decomposing the fractional Fourier transform. Specifically, for each angle a, after some processing steps, only one FFT is needed to obtain the final transform result, which greatly improves the calculation efficiency. On the other hand, in the process of processing the signal, the form of each signal component formed by combination tends to be narrowband or single frequency. From the perspective of computer or processor implementation, the calculation can be realized with higher accuracy, and the calculation accuracy is significantly improved.

[0008] In order to achieve the above-mentioned purpose, the technical solution adopted by the present application is as follows:

[0009] A fractional Fourier transform high-precision fast implementation method for digital signal processing, comprising the following steps:

[0010] S1, according to the actual discretized input signal x s (n) and actual requirements, determining the angle a analysis range of the fractional Fourier transform, the time delay analysis range total size T s of the discretized input signal x m (n), and the time delay analysis step size At m ;

[0011] S2, according to the angle a range, the time delay analysis range total size T m and the time delay analysis step size At m obtained in step S1, calculating the expansion analysis point number N sand the point expansion multiple p, and the discretized input signal x s (n) is zero-padded;

[0012] S3, calculating the digital angular frequency resolution interval and the index value k0 of the start position of the fractional Fourier transform u-axis;

[0013] S4, according to the digital angular frequency resolution interval obtained in step S3, the index value k0, and the discretized input signal x s (n) obtained in step S2, calculating the synthesized signal x' s (n);

[0014] S5, performing anti-aliasing filtering and down-sampling processing on the synthesized signal x' s (n) obtained in step S4, to obtain the down-sampled signal x' sLD (n);

[0015] S6, performing fast Fourier transform (FFT) operation on the down-sampled signal x' sLD (n) obtained in step S5, to obtain the preliminary result X' α (ku s ) of the fractional Fourier transform;

[0016] S7, performing inversion, down-sampling and truncation operations on the preliminary result X' α (ku s ) of the fractional Fourier transform obtained in step S6, to obtain the final result X α (ku s ) of the fractional Fourier transform.

[0017] The step S1 specifically includes:

[0018] For the case that the instantaneous frequency of the discretized input signal x s (n) decreases linearly with time, the analysis range of the angle a should be selected in the interval For the case that the instantaneous frequency of the discretized input signal x s (n) increases linearly with time, the analysis range of the angle a should be selected in the interval ; and the time delay analysis range of the discretized input signal x s (n) is set as [0, T m ], where T m represents the total size of the time delay analysis range, and the time delay analysis step is Δt m , then the total number of time delay analysis points can be obtained as:

[0019]

[0020] Where int(·) represents the rounding operation; for the case where the starting time of the time delay analysis range is not zero, the discretized input signal x is... s (n) Perform advance or delay processing to align the analysis start time to time zero.

[0021] Step S2 specifically involves:

[0022] First, calculate the time delay analysis step size Δt. m The corresponding frequency analysis step size Δf m And the total size T of the time delay analysis range m The corresponding frequency domain analysis length f m Solve using the following formula:

[0023] Δf m =Δt m |cotα| (1-7)

[0024] f m =T m |cotα| (1-8)

[0025] Determine the frequency domain analysis length f m Is it less than or equal to the maximum frequency domain analysis length f? s The sampling rate is used as the reference value. If it is, the process continues; otherwise, it returns to equations (1-7) and (1-8) to reset the total size T of the delay analysis range. m Or angle α;

[0026] The number of extended analysis points N corresponding to the maximum frequency domain range is calculated using the following formula. s :

[0027]

[0028] Where, N s ≥N m N m The point multiplier p represents the total number of points in the time delay analysis compared to the maximum frequency domain range. s Multiples of can be obtained using the following formula:

[0029]

[0030] Among them, symbols This represents rounding up, where N is the discretized input signal x. s (n) The total number of sampling points; after obtaining the point multiplication factor p, the total number of analysis points corresponding to the maximum frequency domain range is pN. s ;

[0031] Then, based on the total number of analysis points pN sFor the discretized input signal x s (n) performs zero-padding, that is, in x s (n) pN is added to the end s -N zeros, so that the total number of points in the signal after zero-padding is exactly pN. s To simplify the description, let the zero-padded discretized input signal still be represented as x. s (n), where the time-domain sequence index n takes values ​​of 0, 1, 2, ..., pN. s -1.

[0032] Step S3 specifically involves:

[0033] Digital angular frequency resolution interval It can be obtained directly from the following formula:

[0034]

[0035] Among them, pN s Total number of analysis points;

[0036] When angle α is in the first quadrant of the time-frequency domain, the index value k0 is calculated using the following formula:

[0037]

[0038] Among them, f bm Indicates the signal to be analyzed, x s The intermediate frequency starting frequency of (n), f bm Satisfying relation: f bm <f s f s Indicates the maximum frequency domain analysis length;

[0039] When angle α is in the fourth quadrant of the time-frequency domain, the index value k0 can be calculated using the following formula:

[0040]

[0041] Where p represents the multiplication factor of the number of points, Δf m This indicates the frequency analysis step size.

[0042] Step S4 specifically involves:

[0043] When angle α is in the first quadrant, x' s (n) is calculated using the following formula:

[0044]

[0045] When angle α is in the fourth quadrant, x' s (n) is calculated using the following formula:

[0046]

[0047] where j represents imaginary unit, n represents time domain sequence index, T s represents sampling interval.

[0048] The step S5 is specifically:

[0049] First, the cutoff angular frequency parameter ω c of the anti-aliasing filter is calculated according to the following formula:

[0050]

[0051] where f m represents frequency domain analysis length, f s represents maximum frequency domain analysis length.

[0052] Then, the down-sampling multiple D is calculated according to the following formula:

[0053]

[0054] where symbol represents rounding down, and max(·) represents taking the maximum value of a group of numbers.

[0055] It is judged whether the down-sampling multiple D is equal to 1. If yes, the anti-aliasing filtering and down-sampling processing are skipped, and the step S6 is directly performed. If no, the anti-aliasing filtering and down-sampling processing are continued.

[0056] Supposing that the length of the anti-aliasing filter is N l , N l is an odd number, the cutoff angular frequency is ω c , the anti-aliasing filter needs to be advancedly shifted to place the symmetry center of the filter at the 0 moment, i.e. at the sequence index n = 0, so as to ensure that the anti-aliasing filter will not produce additional delay to the signal, and the anti-aliasing filter after the above adjustment is denoted as L(n).

[0057] Next, the synthesized signal x' s (n) is subjected to anti-aliasing filtering processing, and the filtered signal x' sL (n) is obtained as shown in the following formula:

[0058]

[0059] where represents circular convolution, and the corresponding circular period is the total analysis point number pN s . The length of the filtered signal x' sL (n) remains unchanged as the total analysis point number pN s , and the time domain sequence index n = 0, 1, 2,..., pN s - 1.

[0060] Then the filtered signal x' sL (n) is down-sampled, and the down-sampled signal x' sLD (n) is expressed as:

[0061] x' sLD (n) = x' sL (nD) (1-17)

[0062] wherein D represents a down-sampling multiple, after the down-sampling processing, the value range of the time domain sequence index n is shortened, and the value becomes

[0063] The step S6 is specifically:

[0064] The preliminary result of the fractional Fourier transform is obtained according to the following formula:

[0065]

[0066] wherein j represents an imaginary unit, u s represents a u-axis discretization step, T s represents a sampling interval, D represents a down-sampling multiple, k' and k are respectively the transform domain sequence indexes before and after correction, k=k'-k0 when the angle α is in the first quadrant, and k=k'+k0 when the angle α is in the fourth quadrant, and the value of k0 is The u-axis discretization step u s is obtained by the following formula:

[0067]

[0068] wherein p represents a point expansion multiple, Δf m represents a frequency analysis step, and it is noted that the signal x' sLD (n) in formula (1-18) is x' sL (n) if the step S5 is not performed.

[0069] The step S7 is specifically:

[0070] Firstly, it is judged according to the angle α whether the preliminary result X' α (ku s ) of the fractional Fourier transform needs to be taken reversely, if the angle α is in the first quadrant, no reversal is needed, and if the angle α is in the fourth quadrant, the preliminary result X' α (ku s ) of the fractional Fourier transform needs to be taken reversely along the u-axis, and the specific operation is performed according to the following formula:

[0071]

[0072] wherein pNs D is the down-sampling factor, and

[0073] X α (ku s ) is then subjected to down-sampling processing, and the down-sampling factor is the point expansion factor p, and the processing is performed according to the following formula:

[0074] X α (ku s ) = X α (kpu s ) (1-21)

[0075] The result X α (ku s ) of formula (1-21) is finally subjected to truncation processing, and the first N m points are extracted, and thus the final fractional Fourier transform result X α (ku s ) is obtained, wherein the modified transform domain sequence index k has a value of k = 0, 1, 2,..., N m -1.

[0076] A high-precision fast implementation system for fractional Fourier transform for digital signal processing, which adopts the processing principle based on the high-precision fast implementation method for fractional Fourier transform for digital signal processing, and comprises:

[0077] A synthetic signal generation module, which determines the analysis range of angle a of the fractional Fourier transform, the total size T m of the time delay analysis range of the signal to be analyzed, and the time delay analysis step size At m according to the actual signal to be processed and actual requirements, and calculates the expansion analysis point number N s corresponding to the maximum frequency domain range and the point expansion factor p according to the angle a range, the total size T m of the time delay analysis range, and the time delay analysis step size At m , and performs zero padding processing on the discretized input signal x s (n), and then calculates the digital angular frequency resolution interval and the index value k0 of the starting position of the fractional Fourier transform u-axis; finally, the synthetic signal x' s (n) is calculated according to the digital angular frequency resolution interval the index value k0, and the discretized input signal x s (n).

[0078] An anti-aliasing filtering and down-sampling module, which performs anti-aliasing filtering and down-sampling processing on the synthetic signal x' s (n) to obtain the down-sampled signal x' sLD (n).

[0079] a fast Fourier transform (FFT) module, for performing fast Fourier transform on the down-sampled signal x' sLD (n) to obtain a preliminary result of fractional Fourier transform X α (ku s );

[0080] a post-processing module for fractional Fourier transform, for performing inverse, down-sampling and truncation operations on the preliminary result of fractional Fourier transform X α (ku s ) to obtain a final result of fractional Fourier transform X α (ku s )。

[0081] An apparatus for high-precision fast implementation of fractional Fourier transform for digital signal processing, comprising:

[0082] a memory for storing a computer program for implementing the method for high-precision fast implementation of fractional Fourier transform for digital signal processing;

[0083] a processor for implementing the method for high-precision fast implementation of fractional Fourier transform for digital signal processing when executing the computer program.

[0084] Compared with the prior art, the present application has the following advantages:

[0085] 1. In the process of calculating the fractional Fourier transform, the present application calculates a composite signal x' s (n) according to the discrete input signal x s (n) to be processed, which realizes narrowband or single frequency processing of the input signal compared with the processing mode of directly using x s (n) for processing, avoids large number calculation in the process of digital calculation for the field of wideband or ultrawideband signal processing, and makes the fractional Fourier transform result more accurate.

[0086] 2. The present application realizes significant reduction of data amount under lossless data information through a series of operations of anti-aliasing filtering, down-sampling and fast Fourier transform (FFT), and has the effect of faster calculation.

[0087] 3. The present application uses a modified transform domain sequence index k as the index of the final fractional Fourier transform result X α (ku s ), so that the starting point index of the fractional Fourier transform corresponds to the starting point of the time delay analysis range of the input signal to be analyzed, has the feature of consistent index correspondence before and after signal transformation, and facilitates actual use.

[0088] 4. The application calculates the extension point number and angular frequency resolution interval of the frequency domain range according to the input signal time delay analysis range and time delay analysis step length of the design requirement, and the design process and result strictly meet the design requirement without additional subsequent processing, thereby increasing the flexibility of the design.

[0089] In summary, the application has the characteristics of high calculation accuracy, fast calculation, and convenient and flexible use. BRIEF DESCRIPTION OF DRAWINGS

[0090] Figure 1 A flow chart of the high-precision fast implementation method of the fractional Fourier transform for digital signal processing according to the application.

[0091] Figure 2 A simulation result graph of the fractional Fourier transform when the search angle is the first one according to the application.

[0092] Figure 3 A simulation result graph of the fractional Fourier transform when the search angle is the 21st one according to the application.

[0093] Figure 4 A simulation result graph of the alpha-u domain two-dimensional distribution of the fractional Fourier transform at all search angles according to the prior art.

[0094] Figure 5 A simulation result graph of the alpha-u domain two-dimensional distribution of the fractional Fourier transform at all search angles according to the application. DETAILED DESCRIPTION

[0095] The technical solutions in the embodiments of the application will be described clearly and completely below with reference to the drawings and specific embodiments. Obviously, the described embodiments are only part of the embodiments of the application, rather than all the embodiments. Based on the embodiments in the application, all other embodiments obtained by those skilled in the art without creative labor fall within the protection scope of the application.

[0096] The technical solutions in the embodiments of the application will be described clearly and completely below with reference to the drawings and specific embodiments. Obviously, the described embodiments are only part of the embodiments of the application, rather than all the embodiments. Based on the embodiments in the application, all other embodiments obtained by those skilled in the art without creative labor fall within the protection scope of the application.

[0097] The definition formula of the fractional Fourier transform (FrFT) to which the application is directed is as follows:

[0098]

[0099] Where x(t) is the signal to be analyzed, t is the time variable, α is the angle, u can be called the displacement factor, and α and u together form the two-dimensional transform domain, K α (t,u) is the FrFT transform kernel, and its expression is:

[0100]

[0101] Where δ(t) represents the unit impulse signal, l is an integer, and j is the imaginary unit. Note that in equation (1-2), the cases α = 2lπ and α = (2l+1)π are only for ensuring the rigor of the fractional Fourier transform expression. In practical applications, the range of values ​​for α is usually set as: That is, α is located in the first or fourth quadrant of the two-dimensional time-frequency domain. The cases where α is located in the second or third quadrant of the time-frequency domain can be easily equivalent to the first or fourth quadrant, so they will not be discussed here.

[0102] The discretized form of equation (1-1) is:

[0103]

[0104] Where n is the index number of the time-domain discrete sequence, and its value is n = 0, 1, 2, ..., N-1, where N is the discretized input signal x. s (n) is the total number of sampling points; u s Let k'u be the discretization step size of the displacement factor u. Then u can be discretized as k'u. s Where k' is the transform domain sequence index number, and its value is related to the actual analysis range of u; transform kernel K α (nT s ,ku s The expression is:

[0105]

[0106] The cases of α = 2lπ and α = (2l+1)π are ignored in the formula; the discretized signal x in formula (1-3) s Let (n) be the sampling form of the signal x(t). Without loss of generality, let x(t) be the intermediate frequency signal after down-conversion and other preprocessing, with its starting time being time 0. Then x s (n) can be expressed as:

[0107]

[0108] In the formula, T s f is the sampling interval. s The sampling rate is the first digit of the sampling rate, and the two are reciprocals of each other.

[0109] like Figure 1A high-precision fast implementation method of fractional Fourier transform for digital signal processing is shown, comprising the following steps:

[0110] S1, according to the actual discrete input signal x s (n) case and actual demand, determine the analysis range of angle α of fractional Fourier transform, the time delay analysis range total size T s (n) of discrete input signal x m , and the time delay analysis step Δt m ;

[0111] S2, according to the angle α range, the time delay analysis range total size T m and the time delay analysis step Δt m obtained in step S1, calculate the expansion analysis point number N s corresponding to the maximum frequency domain range and the point expansion multiple p, and carry out zero padding processing on the discrete input signal x s (n);

[0112] S3, calculate the digital angular frequency resolution interval and the index value k0 of the starting position of fractional Fourier transform u axis;

[0113] S4, according to the digital angular frequency resolution interval obtained in step S3, the index value k0, and the discrete input signal x s (n) obtained in step S2, calculate the synthesis signal x' s (n);

[0114] S5, anti-aliasing filtering and down-sampling processing are carried out on the synthesis signal x' s (n) obtained in step S4, to obtain the down-sampled signal x' sLD (n);

[0115] S6, fast Fourier transform (FFT) operation is carried out on the down-sampled signal x' sLD (n) obtained in step S5, to obtain the preliminary result X' α (ku s ) of fractional Fourier transform;

[0116] S7, reverse, down-sampling and truncation operations are carried out on the preliminary result X' α (ku s ) of fractional Fourier transform obtained in step S6, to obtain the final result X α (ku s ) of fractional Fourier transform.

[0117] The step S1 is specifically:

[0118] for the discrete input signal xs (n) for the case of linearly decreasing instantaneous frequency with time, the analysis range of angle a should be chosen in the interval s (n) for the case of linearly increasing instantaneous frequency with time, the analysis range of angle a should be chosen in the interval ; the time delay analysis range of the discretized input signal x s (n) is set to be [0, T m ], where T m represents the total size of the time delay analysis range, and the time delay analysis step is Δt m , then the total number of time delay analysis points can be obtained as:

[0119]

[0120] where int(·) represents the rounding operation; for the case that the start time of the time delay analysis range is not zero, the discretized input signal x s (n) is processed in advance or delayed to align the analysis start time to the zero time.

[0121] The step S2 is specifically:

[0122] First, the time delay analysis step Δt m , the corresponding frequency analysis step Δf m , and the total size of the time delay analysis range T m , the corresponding frequency domain analysis length f m are solved according to the following formula:

[0123] Δf m = Δt m | cot a | (1-7)

[0124] f m = T m | cot a | (1-8)

[0125] It is judged whether the frequency domain analysis length f m is less than or equal to the maximum frequency domain analysis length f s , i.e. the sampling rate, if yes, continue to execute, if no, return to formula (1-7) and formula (1-8) to reset the total size of the time delay analysis range T m or the angle a;

[0126] The extended analysis point number N s corresponding to the maximum frequency domain range is calculated according to the following formula:

[0127]

[0128] where N s ≥ N​m , N m is the total number of points of the input signal, and p is the point expansion factor, which represents the ratio of the total number of points N of the input signal to the number of analysis points N of the maximum frequency range s , which is calculated as follows:

[0129]

[0130] wherein the symbol represents rounding up, and N is the total number of sampling points of the discretized input signal x s (n); after obtaining the point expansion factor p, the total number of analysis points corresponding to the maximum frequency range is pN s .

[0131] Then, the discretized input signal x s (n) is zero-padded according to the total number of analysis points pN s , i.e. pN s -N zeros are appended to the end of x s (n) so that the total number of points of the zero-padded signal is exactly pN s , for simplicity of representation, the zero-padded discretized input signal is still denoted as x s (n), wherein the time-domain sequence index n takes values of 0, 1, 2,..., pN s -1.

[0132] The step S3 is specifically:

[0133] The digital angular frequency resolution interval is directly obtained by the following formula:

[0134]

[0135] wherein pN s is the total number of analysis points;

[0136] For the calculation of the index value k0of the starting position of the fractional Fourier transform u-axis, it is necessary to discuss according to the quadrant of the time-frequency domain where the angle a is located:

[0137] When the angle a is located in the first quadrant of the time-frequency domain, the index value k0is calculated as follows:

[0138]

[0139] wherein f bm represents the intermediate frequency starting frequency of the signal to be analyzed x s (n), and f bm satisfies the relationship: f bm <f s , f s represents the maximum frequency domain analysis length;

[0140] When the angle a is in the fourth quadrant of the time-frequency domain, the index value k0 can be calculated as follows:

[0141]

[0142] where p represents a point expansion multiple, Δf represents a frequency analysis step size, and f represents a frequency domain analysis length. m

[0143] The step S4 is specifically:

[0144] When the angle a is in the first quadrant, x' is calculated as follows: s (n) is calculated as follows:

[0145]

[0146] When the angle a is in the fourth quadrant, x' is calculated as follows: s (n) is calculated as follows:

[0147]

[0148] where j represents an imaginary unit, n represents a time domain sequence index, T represents a sampling interval, and f represents a frequency domain analysis length. s

[0149] The step S5 is specifically:

[0150] First, a cutoff angular frequency parameter ω of the anti-aliasing filter is calculated according to the following formula: c

[0151]

[0152] where f represents a frequency domain analysis length, f represents a maximum frequency domain analysis length, and f represents a sampling interval. m s

[0153] Then, a down-sampling multiple D is calculated according to the following formula:

[0154]

[0155] where the symbol represents a floor function, and max(·) represents a maximum value of a set of numbers.

[0156] It is determined whether the down-sampling multiple D is equal to 1. If yes, the anti-aliasing filtering and down-sampling processing below are skipped, and the step S6 is directly performed. If no, the anti-aliasing filtering and down-sampling processing below is continued to be performed.

[0157] The anti-aliasing filter should be a linear phase FIR low-pass filter designed according to a known technique, and the filter length should be set as an odd number, and the cutoff angular frequency of the filter is set as ω in the formula (1-14).​​​​​​c In this embodiment, a linear phase window function low-pass filter is used as an anti-aliasing filter.

[0158] Suppose the length of the designed anti-aliasing filter is N l , N l is an odd number, and the cutoff angular frequency is ω c , the anti-aliasing filter needs to be advancedly shifted so that the symmetric center of the filter is at time 0, i.e., at sequence index n = 0, to ensure that the anti-aliasing filter does not cause additional delay to the signal. Suppose the anti-aliasing filter after the above adjustment is L(n).

[0159] Next, the synthesized signal x' s (n) is subjected to anti-aliasing filter processing to obtain the filtered signal x' sL (n) as shown in the following formula:

[0160]

[0161] wherein, represents circular convolution, and the corresponding circular period is the total number of analysis points pN s The length of the filtered signal x' sL (n) remains unchanged at the total number of analysis points pN s , and the time domain sequence index n = 0, 1, 2,..., pN s -1.

[0162] Then, the filtered signal x' sL (n) is subjected to downsampling processing, and the downsampled signal x' sLD (n) is represented as:

[0163] x' sLD (n) = x' sL (nD) (1-17)

[0164] wherein, D represents the downsampling multiple, and after the downsampling processing, the value range of the time domain sequence index n is shortened, and the value becomes

[0165] The step S6 is specifically:

[0166] The preliminary result of the fractional Fourier transform is obtained according to the following formula:

[0167]

[0168] wherein, j represents the imaginary unit, u s represents the u-axis discretization step, and T srepresents a sampling interval, D represents a down-sampling multiple, k' and k are respectively a transform domain sequence index before and after correction, k=k'-k0 when the angle a is in the first quadrant, k=k'+k0 when the angle a is in the fourth quadrant, and the value of k0 is The u-axis discretization step u in the formula s is obtained by the following formula:

[0169]

[0170] Wherein, p represents a point expansion multiple, Δf m represents a frequency analysis step, and it is noted that if the step S5 is not performed (corresponding to the case of D=1), the signal x' sLD (n) in the formula (1-18) is x' sL (n).

[0171] The step S7 is specifically:

[0172] First, according to the angle a, it is judged whether the preliminary result X' α (ku s ) of the fractional Fourier transform needs to be reversed, if the angle a is in the first quadrant, no reversal is needed, if the angle a is in the fourth quadrant, a reversal operation is needed to be performed to the preliminary result X' α (ku s ) of the fractional Fourier transform along the u-axis, which is specifically performed according to the following formula:

[0173]

[0174] Wherein, pN s is a total analysis point number, and D is a down-sampling multiple;

[0175] Then, the down-sampling processing is continuously performed to X" α (ku s ), and the down-sampling multiple is the point expansion multiple p, which is performed according to the following formula:

[0176] X α (ku s ) = X" α (kpu s ) (1-21)

[0177] Finally, the truncation processing is performed to the result X α (ku s ) of the formula (1-21), and the first N m points are extracted, and thus, the final fractional Fourier transform result X α (ku s ) is obtained, wherein the transform domain sequence index k after correction is k=0, 1, 2,..., N m -1.

[0178] In principle, the present application makes a new form of decomposition on the fractional Fourier transform expression, and obtains a new decomposition form different from the prior art:

[0179]

[0180] From the expression form, the superposition and part of it Only one fast Fourier transform (FFT) operation can obtain the fractional Fourier transform result of a certain signal under an angle α, without the need for multiple fast Fourier transform (FFT) and inverse fast Fourier transform (IFFT) operations as in the prior art. However, if the FFT operation is directly performed on the superposition and part, the number of points is pN s , which is too large and not efficient. Therefore, in the present application, the synthesized signal x' s (n) is further subjected to anti-aliasing filtering and down-sampling processing to further reduce the number of signal points. Then, the FFT operation is performed on the processed synthesized signal, which greatly improves the calculation efficiency.

[0181] On the other hand, it can be analyzed that the synthesized signal x' s (n) in the present application (its form is shown in (1-13a) and (1-13b)) realizes narrowband or single frequency with respect to the original signal x s (n), which can realize higher precision calculation in actual processing. As a comparison, the prior art usually involves mutual operation between two or more synthesized signals or constructed signals in the form of wideband linear frequency modulation, which will inevitably produce errors in actual numerical calculation process, and the wider the frequency band of the signal to be processed and the longer the time length, the greater the processing error. This reflects the advantage of high calculation precision of the present application in the field of wideband or ultrawideband signal processing.

[0182] A fractional Fourier transform high-precision fast implementation system for digital signal processing, whose processing principle adopts the fractional Fourier transform high-precision fast implementation method for digital signal processing, comprising:

[0183] A synthesized signal generation module determines the analysis range of the angle α of the fractional Fourier transform, the total size T m of the time delay analysis range of the signal to be analyzed, and the time delay analysis step size Δt m according to the actual signal to be processed and the actual demand, and then calculates the expansion analysis point number N s corresponding to the maximum frequency domain range and the point expansion multiple p according to the angle α range, the total size T m of the time delay analysis range, and the time delay analysis step size Δt m , and performs the expansion analysis on the discretized input signal x s(n) zero padding and then calculating the digital angular frequency resolution interval and the index value k0 of the start position of the fractional Fourier transform u-axis; finally, according to the digital angular frequency resolution interval the index value k0 and the discretized input signal x s (n) calculating the synthesis signal x' s (n) for implementing steps S1-S4 of the high-precision fast implementation method of the fractional Fourier transform for digital signal processing;

[0184] the anti-aliasing filtering and down-sampling module performs anti-aliasing filtering and down-sampling processing on the synthesis signal x' s (n) to obtain the down-sampled signal x' sLD (n) for implementing step S5 of the high-precision fast implementation method of the fractional Fourier transform for digital signal processing;

[0185] the fast Fourier transform (FFT) module performs fast Fourier transform (FFT) operation on the down-sampled signal x' sLD (n) to obtain the preliminary result X' α (ku s ) of the fractional Fourier transform;

[0186] the fractional Fourier transform post-processing module performs inversion, down-sampling and truncation operation on the preliminary result X' α (ku s ) of the fractional Fourier transform to obtain the final result X α (ku s ) of the fractional Fourier transform;

[0187] A high-precision fast implementation device of the fractional Fourier transform for digital signal processing, comprising:

[0188] a memory for storing a computer program for implementing the high-precision fast implementation method of the fractional Fourier transform for digital signal processing according to any one of claims 1-8;

[0189] a processor for implementing the high-precision fast implementation method of the fractional Fourier transform for digital signal processing according to any one of claims 1-8 when executing the computer program.

[0190] The present application is directed to a high-precision fast calculation method of the fractional Fourier transform, which is most commonly used in the field of signal processing or system, and particularly relates to specific fields such as radar signal processing, radio frequency signal processing and digital signal processing.

[0191] The above merely describes preferred embodiments of the present application but should not be used to limit the protective scope of the present application. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present application shall fall within the protective scope of the present application.

[0192] The effects of the present application can be further illustrated by the following simulation:

[0193] 1. Simulation conditions

[0194] Suppose that the signal to be analyzed is a wideband linear frequency modulation pulse signal, the signal pulse width is 0.5ms, the bandwidth is 500MHz, the modulation mode is linear increase of frequency with time, the angle a corresponding to the time-frequency domain is in the fourth quadrant, the input signal-to-noise ratio is-20dB, the intermediate frequency sampling rate is set to 1.2GHz, the total sampling points of the signal to be analyzed are 8x10 5 , the total size of the time delay analysis is 66.7us, the time delay analysis step is 1.5ns, the total number of search angles is set to 60, and the real angle a corresponding to the signal to be analyzed is in the 21st search angle sequence. The simulation computer configuration: CPU 3.8GHz / 24 cores, memory 16GB. The parameters here have no substantial influence on the simulation results.

[0195] 2. Simulation content

[0196] Under the above simulation conditions, the fractional Fourier transform processing is performed on the signal to be analyzed by using the fractional Fourier transform high-precision fast implementation method for digital signal processing according to the present application and the prior art respectively, and the fractional Fourier transform results are recorded as shown in Figures 2-5 , and the simulation running time consumption results are shown in Table 1.

[0197] Table 1. Comparison of operation time consumption of prior art and method of the present application

[0198] Category Prior art Method of the invention Calculation time (seconds) 52.37 12.28

[0199] Figure 2 And Figure 3 The fractional Fourier transform results under the conditions of incomplete angle matching and complete angle matching are given respectively, and from the results, it can be seen that the peak energy corresponding to the fractional Fourier transform high-precision fast implementation method for digital signal processing according to the present application is more concentrated, and the amplitude is also higher, while the peak energy corresponding to the prior art is more dispersed, and the amplitude is slightly lower, which shows that the calculation precision of the fractional Fourier transform high-precision fast implementation method for digital signal processing according to the present application is higher, and a higher matching correlation peak can be obtained. Figure 4 And Figure 5The fractional Fourier transform two-dimensional distribution result of the score further embodies the peak energy divergence effect caused by the loss of calculation accuracy of the prior art, and the high-precision fast implementation method for fractional Fourier transform for digital signal processing realizes high-precision calculation and obtains a more concentrated correlation peak.

[0200] Table 1 gives the total time consumption of the fractional Fourier transform under 60 search angles, and it can be obviously seen that the time consumption of the high-precision fast implementation method for fractional Fourier transform for digital signal processing is 12.28s, which is much smaller than 52.37s of the prior art, embodying the high efficiency of the high-precision fast implementation method for fractional Fourier transform for digital signal processing in realizing fast calculation of the fractional Fourier transform.

[0201] It should be noted that the signal to be analyzed in the simulation example belongs to a high-bandwidth and high-pulse-width signal, and the prior art has obvious loss of calculation accuracy and high time consumption, which is the main application scenario in which the advantages of the high-precision fast implementation method for fractional Fourier transform for digital signal processing are embodied. For other application scenarios, the quantitative conclusions of the simulation example are not applicable.

Claims

1. A high-precision and fast method for implementing fractional Fourier transform in digital signal processing, characterized in that: Includes the following steps: S1, based on the actual discretized input signal x s (n) Based on the situation and actual needs, determine the analytical range of the angle α of the fractional Fourier transform and the discretized input signal x. s Total size T of the time delay analysis range of (n) m and the time delay analysis step size Δt m ; S2, based on the angle α range and the total size T of the time delay analysis range obtained in step S1. m Time delay analysis step size Δt m Calculate the number of extended analysis points N corresponding to the maximum frequency domain range. s And the point expansion factor p, and the discretized input signal x s (n) is padded with zeros; S3, Calculate the digital angular frequency resolution interval And the index value k0 of the starting position of the fractional Fourier transform u-axis; S4, Based on the digital angular frequency resolution interval obtained in step S3 The index value k0, and the discretized input signal x obtained in step S2 s (n) Calculate the synthesized signal x' s (n); S5, the synthesized signal x' obtained in step S4 s (n) Perform anti-aliasing filtering and downsampling to obtain the downsampled signal x' sLD (n); S6, the downsampled signal x' obtained in step S5 sLD (n) Perform a Fast Fourier Transform (FFT) operation to obtain the preliminary result X' of the Fractional Fourier Transform. α (ku s ); S7, the preliminary result X' of the fractional Fourier transform obtained in step S6. α (ku s Perform inversion, downsampling, and truncation operations to obtain the final fractional Fourier transform result X. α (ku s ).

2. The method for high-precision and fast implementation of fractional Fourier transform for digital signal processing according to claim 1, characterized in that: Step S1 specifically involves: For discretized input signal x s (n) represents the case where the instantaneous frequency decreases linearly with time; the analysis range for angle α should be within the interval [interval missing]. Selected from the discretized input signal x s (n) represents the case where the instantaneous frequency increases linearly with time, and the analysis range of angle α should be within the interval [missing information]. Select from the options; set the discretized input signal x. s The time delay analysis range of (n) is [0, T m ], where T m This represents the total size of the time delay analysis range, with a time delay analysis step size of Δt. m Then the total number of delay analysis points can be obtained as follows: Where int(·) represents the rounding operation; for the case where the starting time of the time delay analysis range is not zero, the discretized input signal x is... s (n) Perform advance or delay processing to align the analysis start time to time zero.

3. The method for high-precision and fast implementation of fractional Fourier transform for digital signal processing according to claim 1, characterized in that: Step S2 specifically involves: First, calculate the time delay analysis step size Δt. m The corresponding frequency analysis step size Δf m And the total size T of the time delay analysis range m The corresponding frequency domain analysis length f m Solve using the following formula: Δf m =Δt m |cotα| (1-7) f m =T m |cotα| (1-8) Determine the frequency domain analysis length f m Is it less than or equal to the maximum frequency domain analysis length f? s The sampling rate is used as the reference value. If it is, the process continues; otherwise, it returns to equations (1-7) and (1-8) to reset the total size T of the delay analysis range. m Or angle α; The number of extended analysis points N corresponding to the maximum frequency domain range is calculated using the following formula. s : Where, N s ≥N m N m The point multiplier p represents the total number of points in the time delay analysis compared to the maximum frequency domain range. s Multiples of can be obtained using the following formula: Among them, symbols This represents rounding up, where N is the discretized input signal x. s (n) The total number of sampling points; after obtaining the point multiplication factor p, the total number of analysis points corresponding to the maximum frequency domain range is pN. s ; Then, based on the total number of analysis points pN s For the discretized input signal x s (n) performs zero-padding, that is, in x s (n) pN is added to the end s -N zeros, so that the total number of points in the signal after zero-padding is exactly pN. s To simplify the description, let the zero-padded discretized input signal still be represented as x. s (n), where the time-domain sequence index n takes values ​​of 0, 1, 2, ..., pN. s -1.

4. The method for high-precision and fast implementation of fractional Fourier transform for digital signal processing according to claim 1, characterized in that: Step S3 specifically involves: Digital angular frequency resolution interval It can be obtained directly from the following formula: Among them, pN s Total number of analysis points; When angle α is in the first quadrant of the time-frequency domain, the index value k0 is calculated using the following formula: Among them, f bm Indicates the signal to be analyzed, x s The intermediate frequency starting frequency of (n), f bm Satisfying relation: f bm <f s f s Indicates the maximum frequency domain analysis length; When angle α is in the fourth quadrant of the time-frequency domain, the index value k0 can be calculated using the following formula: Where p represents the multiplication factor of the number of points, Δf m This indicates the frequency analysis step size.

5. A high-precision and fast implementation method for fractional Fourier transform in digital signal processing according to claim 1, characterized in that: Step S4 specifically involves: When angle α is in the first quadrant, x' s (n) is calculated using the following formula: When angle α is in the fourth quadrant, x' s (n) is calculated using the following formula: Where j represents the imaginary unit, n represents the time-domain sequence index, and T s Indicates the sampling interval.

6. The method for high-precision and fast implementation of fractional Fourier transform for digital signal processing according to claim 1, characterized in that: Step S5 specifically involves: First, calculate the cutoff angular frequency parameter ω of the anti-aliasing filter according to the following formula. c : Among them, f m f represents the frequency domain analysis length. s Indicates the maximum frequency domain analysis length; Then calculate the downsampling factor D according to the following formula: Among them, symbols This indicates rounding down, while max(·) indicates taking the maximum value of a set of numbers; Determine if the downsampling factor D is equal to 1. If it is, skip the anti-aliasing filtering and downsampling process and proceed directly to step S6. If not, continue to execute the anti-aliasing filtering and downsampling process. Let the length of the anti-aliasing filter be N. l N l It is an odd number, and the cutoff angular frequency is ω. c Therefore, the anti-aliasing filter needs to be shifted ahead so that the center of symmetry of the filter is placed at time 0, that is, at the corresponding sequence index n=0, to ensure that the anti-aliasing filter will not cause additional delay to the signal. Let the anti-aliasing filter adjusted as above be L(n). Next, the synthesized signal x' s (n) Perform anti-aliasing filtering to obtain the filtered signal x' sL (n) is shown in the following formula: in, This represents a circular convolution, with a corresponding cycle period of pN (the total number of analysis points). s Filtered signal x' sL The length of (n) remains constant as the total number of analysis points pN. s The time-domain sequence index remains unchanged, n = 0, 1, 2, ..., pN s -1; Then the filtered signal x' sL (n) Perform downsampling processing on the signal x' after downsampling. sLD (n) is represented as: x' sLD (n)=x' sL (nD) (1-17) Where D represents the downsampling factor. After downsampling, the range of values ​​for the time-domain sequence index n is shortened, and its value becomes...

7. A high-precision and fast implementation method for fractional Fourier transform in digital signal processing according to claim 1, characterized in that: Step S6 specifically involves: The preliminary results of the fractional Fourier transform are obtained according to the following formula: Where j represents the imaginary unit, u s T represents the step size for discretization along the u-axis. s The sampling interval is represented by D, the downsampling factor is represented by k' and k', and the transform domain sequence indices before and after correction are respectively. When the angle α is in the first quadrant, k = k' - k0, and when the angle α is in the fourth quadrant, k = k' + k0. The values ​​are... The discretization step size u of the u-axis in the formula s We obtain it from the following formula: Where p represents the multiplication factor of the number of points, Δf m This represents the frequency analysis step size. Note that if step S5 is not performed, then the signal x' in equation (1-18) will be... sLD (n) is x' sL (n).

8. A high-precision and fast implementation method for fractional Fourier transform in digital signal processing according to claim 1, characterized in that: Step S7 specifically involves: First, determine the preliminary result X' of the fractional Fourier transform based on angle α. α (ku s Whether or not inversion is needed depends on the angle α. If α is in the first quadrant, inversion is not needed. If α is in the fourth quadrant, inversion is needed for the preliminary result X' of the fractional Fourier transform. α (ku s Perform the inversion operation along the u-axis, specifically as follows: Among them, pN s Where D is the total number of analysis points, and D is the downsampling factor. Then for X” α (ku s Continue with downsampling. The downsampling factor is the point multiplication factor p, calculated using the following formula: X α (stand s )=X” α (cap s ) (1-21) Finally, the result X of equation (1-21) α (ku s Perform truncation and extract the first N. m At this point, the final fractional Fourier transform result X is obtained. α (ku s The corrected transform domain sequence index k takes the values: k = 0, 1, 2, ..., N m -1.

9. A high-precision, fast implementation system for fractional Fourier transform in digital signal processing, wherein the processing principle adopts the high-precision, fast implementation method for fractional Fourier transform in digital signal processing as described in any one of claims 1 to 8, characterized in that: include: The synthesized signal generation module determines the analysis range of the angle α of the fractional Fourier transform and the total size T of the time delay analysis range of the signal to be analyzed, based on the actual signal conditions and actual requirements. m and the time delay analysis step size Δt m Then, based on the range of angle α and the total size T of the time delay analysis range, m Time delay analysis step size Δt m Calculate the number of extended analysis points N corresponding to the maximum frequency domain range. s And the point expansion factor p, and the discretization of the input signal x. s (n) Perform zero-padding, and then calculate the digital angular frequency resolution interval. and the index value k0 of the starting position of the fractional Fourier transform u-axis; finally, the interval is resolved based on the digital angular frequency. Index value k0, and discretized input signal x s (n) Calculate the synthesized signal x' s (n); The anti-aliasing filtering and downsampling module processes the synthesized signal x' s (n) Perform anti-aliasing filtering and downsampling to obtain the downsampled signal x' sLD (n); The Fast Fourier Transform (FFT) module performs a fast Fourier transform on the downsampled signal x'. sLD (n) Perform a Fast Fourier Transform (FFT) operation to obtain the preliminary result X' of the Fractional Fourier Transform. α (ku s ); The fractional Fourier transform post-processing module processes the initial result X' of the fractional Fourier transform. α (ku s Perform inversion, downsampling, and truncation operations to obtain the final fractional Fourier transform result X. α (ku s ).

10. A device for high-precision and rapid implementation of fractional Fourier transform for digital signal processing, characterized in that, include: Memory: Used to store a computer program that implements the high-precision and fast implementation method of fractional Fourier transform for digital signal processing as described in any one of claims 1 to 8; Processor: Used to implement, when executing the computer program, a high-precision and fast method for implementing fractional Fourier transform for digital signal processing as described in any one of claims 1 to 8.

Images