CN101303689A - Method for implementing fractional order Fourier transformation based on multi-sample - Google Patents

Abstract

The invention relates to a fractional Fourier transform realizing method based on multi-sampling, belongs to the field of signal processing and can be used for improving the calculation efficiency realized in real time. The method decomposes the original arithmetic into two parallel multi-phase equivalent structures by adopting a multi-sampling signal processing principle based on a discrete sampling arithmetic, thus removing redundant operations, optimizing the arithmetic flows and solves the problems of slow speed, large time lag and not being not beneficial to real time processing of an Ozaktas discrete sampling arithmetic to realize real time calculation; compared with the original arithmetic, the flow calculation load is smaller and the efficiency is higher; besides, the parallel structures are more suitable for the appearance of hardware; moreover, an more effective window-adding semi-band filter is selected at the part of inputting the interpolation filter of the pre-treatment of a sequence, thus not only solving the edge oscillation problem occurred during low orders, but also being beneficial to pipeline processing and improving the calculation efficiency when realizing as the length of the filter is shorter.

Description

A kind of fraction Fourier conversion implementation method based on multisampling

Technical field

The invention belongs to the signal Processing field, be specifically related to the efficient implementation method of a kind of fraction Fourier conversion, improve the counting yield of real-time implementation based on multisampling.

Background technology

Fraction Fourier conversion has widespread use at optical field at first, and Almeida was interpreted as the rotation of signal at time-frequency plane to fraction Fourier conversion in 1993, was the popularization of classical Fourier transform; After Ozaktas had proposed a kind of discrete sampling type algorithm suitable with the FFT computing velocity in 1996, fraction Fourier conversion just began to be applied in the signal Processing field.Fraction Fourier conversion can be regarded a kind of unified time-frequency conversion as, reflected simultaneously signal the time, the information of frequency domain, different with quadratic form time-frequency distributions commonly used is that it represents time-frequency information with unitary variant, and there is not the cross term puzzlement, compare with traditional Fourier transform (being a special case of fraction Fourier conversion in fact), it is suitable for handling non-stationary signal, especially chirp class signal, and (the conversion exponent number a) for many free parameters, therefore fraction Fourier conversion often can access traditional time-frequency distributions or effect that Fourier transform can not get under certain conditions, and, therefore when obtaining better effect, do not need to pay too many calculation cost because it has the fast discrete algorithm of comparative maturity.There are following 4 kinds of modes in the signal Processing field to the application of fraction Fourier conversion at present:

(1) fraction Fourier conversion is a kind of unified time-frequency conversion, along with exponent number rises to 1 continuously from 0, fraction Fourier conversion shows signal and progressively changes to all changes feature of frequency domain from time domain, and bigger choice can be provided for the time frequency analysis of signal; The most directly utilize that mode is exactly during with tradition, the application of frequency domain to fractional order Fourier domain obtaining the improvement on some performance, as the filtering of fractional order Fourier domain etc.

(2) fraction Fourier conversion can be understood as the decomposition of chirp base, therefore, and its very suitable processing chirp class signal, and chirp class signal often runs at radar, communication, sonar and occurring in nature.Wave beam formation, Target Recognition and frequency sweep that fraction Fourier conversion can be handled in these fields are disturbed inhibition or the like problem.

(3) fraction Fourier conversion is the rotation to time-frequency plane, utilize this point can set up the relation of fraction Fourier conversion and time frequency analyzing tool, both can be used for estimating instantaneous frequency, recover phase information, can be used for designing new time frequency analyzing tool again, be signal spreading method of basis function etc. as: TTFT, with the fraction Fourier conversion of Gaussian function.

(4) compare with Fourier transform, therefore the many free parameters of fraction Fourier conversion can access better effect in some application scenario, as: digital watermarking and image encryption.

The Fourier Transform of Fractional Order definition is as follows:

$F^a\left[f\left(u\right)\right]\≡\left\{F^{a}f\right\}\left(u\right)\≡{\∫}_{-\∞}^{\∞}{B}_a(u,x)f\left(x\right)\mathrm{dx},$

B _a(u，x)≡A _φexp[iπ(u ²cotφ-2ux cscφ+x ² cotφ)]，

A wherein _φ≡ { exp (i π sgn (sin φ)/4+i φ/2) }/| sin φ |, φ ≡ a pi/2.

Rearranging following formula can obtain:

$\left\{F^{a}f\right\}\left(u\right)=A_{\mathrm{\φ}}{e}^{\mathrm{j\π}{u}^2\cot \mathrm{\φ}}{\∫}_{-\∞}^{\∞}{e}^{-j2\mathrm{\π}\mathrm{ux}\csc \mathrm{\φ}}\left[e^{\mathrm{j\π}{x}^2\cot \mathrm{\φ}}f\left(x\right)\right]\mathrm{dx}$

The discrete logarithm of fraction Fourier conversion has a variety of, wherein the fast algorithm of Turk Ozaktas proposition is a kind of of present main flow algorithm, its ultimate principle is that the fraction Fourier conversion through the later signal f (x) of dimension normalization is resolved into following 3 easy steps, and each step is carried out discretize:

1. at first use chirp signal modulation signal f (x) to obtain g (x):

$g\left(x\right)=f\left(x\right)\exp (-\mathrm{j\π}{x}^2\tan \frac{\mathrm{\φ}}{2}).$

2. modulation signal g (x) and another chirp signal convolution:

$G^a\left(u\right)=A_{\mathrm{\φ}}{\∫}_{-\∞}^{\∞}\exp \left[\mathrm{j\π}{(u-x)}^2\csc \mathrm{\φ}\right]g\left(x\right)\mathrm{dx}$

3. with the signal after the chirp signal modulation convolution:

$\left\{F^{a}f\right\}\left(u\right)=G^a\left(u\right)\exp (-\mathrm{j\π}{u}^2\tan \frac{\mathrm{\φ}}{2})$

Can list of references " the dimension Study on Normalization in the fraction Fourier conversion numerical evaluation " about dimension normalization, be published in " Beijing Institute of Technology's journal " fourth phase in 2005, the author is Zhao Xinghao, Deng Bing, happy and carefree).

The specific implementation process of the fast algorithm that Ozaktas proposes is:

It is the center that the Wigner-Ville distribution of putative signal f (x) is limited to initial point, and diameter is in the circle of Δ.Then when fractional order conversion order a satisfy 0.5≤| during a|≤1.5,

Maximum bandwidth be limited in the Δ, be sampling interval with 1/2 Δ, rebuild with shannon formula

Can get:

$e^{\mathrm{j\π}{x}^2\cot \mathrm{\φ}}f\left(x\right)=\underset{n=-N}{\overset{N}{\mathrm{\Σ}}}{e}^{\mathrm{j\π}{\left(\frac{n}{2\mathrm{\Δ}}\right)}^2\cot \mathrm{\φ}}f\left(\frac{n}{2\mathrm{\Δ}}\right)\sin c\left(2\mathrm{\Δ}(x-\frac{n}{2\mathrm{\Δ}})\right)---\left(2\right)$

N=Δ wherein ²(2) formula substitution (1) formula abbreviation is obtained:

$\left\{F^{a}f\right\}\left(u\right)=\frac{A_{\mathrm{\φ}}}{2\mathrm{\Δ}}{e}^{\mathrm{j\π}{u}^2\cot \mathrm{\φ}}\underset{n=-N}{\overset{N}{\mathrm{\Σ}}}{e}^{-j2\mathrm{\πu}\frac{n}{2\mathrm{\Δ}}\csc \mathrm{\φ}}{e}^{\mathrm{j\π}{\left(\frac{n}{2\mathrm{\Δ}}\right)}^2\cot \mathrm{\φ}}f\left(\frac{n}{2\mathrm{\Δ}}\right)---\left(3\right)$

Be sampling interval with 1/2 Δ then, in [Δ/2, Δ/2] scope, discretize carried out in the fractional order territory, order $u=\frac{m}{2\mathrm{\Δ}},$ Then:

$\left\{F^{a}f\right\}\left(\frac{m}{2\mathrm{\Δ}}\right)=\frac{A_{\mathrm{\φ}}}{2\mathrm{\Δ}}{e}^{\mathrm{j\π}{\left(\frac{m}{2\mathrm{\Δ}}\right)}^2\cot \mathrm{\φ}}\underset{n=-N}{\overset{N}{\mathrm{\Σ}}}{e}^{-j2\mathrm{\π}\frac{\mathrm{mn}}{{\left(2\mathrm{\Δ}\right)}^2}\csc \mathrm{\φ}}{e}^{\mathrm{j\π}{\left(\frac{n}{2\mathrm{\Δ}}\right)}^2\cot \mathrm{\φ}}f\left(\frac{n}{2\mathrm{\Δ}}\right)---\left(4\right)$

Wherein-and N≤m≤N, further put following formula in order and obtain:

$\left\{F^{a}f\right\}\left(\frac{m}{2\mathrm{\Δ}}\right)=\frac{A_{\mathrm{\φ}}}{2\mathrm{\Δ}}{e}^{-\mathrm{j\π}{\left(\frac{m}{2\mathrm{\Δ}}\right)}^2\tan \frac{\mathrm{\φ}}{2}}\underset{n=-N}{\overset{N}{\mathrm{\Σ}}}{e}^{j\mathrm{\π}{\left(\frac{m-n}{2\mathrm{\Δ}}\right)}^2\csc \mathrm{\φ}}{e}^{-\mathrm{j\π}{\left(\frac{n}{2\mathrm{\Δ}}\right)}^2\tan \frac{\mathrm{\φ}}{2}}f\left(\frac{n}{2\mathrm{\Δ}}\right)]---\left(5\right)$

Suppose that original signal sequence length is N, the normalization sampling rate $\mathrm{\Δ}=\sqrt{N},$ Sampling interval is 1/ Δ, and g (x) in above 3 steps, G ^a(u) sampling interval all is 1/2 Δ, therefore need the original signal sample sequence is carried out 2 times of interpolations before the first step chirp modulation carrying out, specific practice is earlier with burst null value interpolation, level and smooth by a low-pass filter h (n) then, result and 2N length chirp sequence through the null value interpolation multiply each other, make convolution with the 4N length sequences again, convolution results intercepting 2N point multiplies each other with first chirp sequence and a constant coefficient again, also will carry out 2 times at last and extract so that obtain the transformation results sequence { F that sampling interval is 1/ Δ ^aF (m/ Δ), similar FFT, total conversion process are the N inputs, N output.

There are following two shortcomings in classical implementation method according to above-mentioned algorithm flow:

(1) calculating process speed is slow, and time delay is big.A N input, the conversion process intermediate demand of N output is through the convolution of a 4N length, the FFT of 3 4N length or IFFT conversion just.

(2) use more null value interpolating method at present and mainly contain following two kinds, each have different shortcomings: a kind of sinc of being sequence of function is made the level and smooth null value interpolation of low-pass filter sequence, general this filter length is 4 times of list entries length, effective, but calculated amount is big, be unfavorable for that hardware realizes, and order is when low, vibration appears in the edge of transformation results; Another kind is to use lagrange-interpolation, and this method computing velocity is fast, but effect is bad, and especially the marginal error in result of calculation is big.A.Bultheel compares these two kinds of methods, specifically please see document " Computation of the Fractional Fourier transform ", be published in " Applied and Computational Harmonic Analysis " the 16th phase in 2004.

Summary of the invention

The objective of the invention is in order to solve in the above-mentioned classical implementation method computing velocity slow, time delay is big, be unfavorable for problems such as processing in real time, a kind of fraction Fourier conversion implementation method based on multisampling has been proposed, the relative primal algorithm flow process of this method calculated amount is littler, efficient is higher, and owing to be parallel organization, be more suitable for hardware and realize; In addition, the present invention has partly selected more effective windowing half-band filter for use at the pretreated filtering interpolation of list entries, the edge oscillation problem that has occurred when not only having solved low order, and because filter length is shorter, help stream treatment, improve the counting yield when realizing.

The ultimate principle of the inventive method is on the fast algorithm basis that Ozaktas proposes, adopt the multisampling signal processing theory, primal algorithm is resolved into the parallel heterogeneous equivalent structure of two-way, removed redundant operation, optimize algorithm flow, improved counting yield.

The objective of the invention is to be achieved through the following technical solutions.

A kind of fraction Fourier conversion implementation method based on multisampling of the present invention comprises following 5 steps, and each step comprises the branch road of two concurrent operations:

1. with the heterogeneous component r of input signal sequence f (n/ Δ) by interpolation filter h (n) _i(k) carry out filtering, obtain sequence f _i(m/ Δ), i=0 wherein, 1 an expression way.|n|≤(N-1)/2， $\mathrm{\Δ}=\sqrt{N}$ The normalized sampling interval of expression input signal sequence, N represents input signal sequence length, is odd number, wherein:

H (n)=2w (n) sin (π n/2)/(π n), | n|≤L, 2L+1 are filter length;

w(n)＝0.5·[1+cos(nπ/L)]，|n|≤L；

r ₀(n)＝h(2n+1)，r ₁(n)＝h(2n)＝δ(n)；

2. use two heterogeneous component l of chirp signal chirpA (n/ (2 Δ)) _i(k/ Δ) respectively with the 1. filtering f as a result of respective branch of step _i(m/ Δ) multiplies each other and obtains g _i(k/ Δ), wherein:

$\mathrm{chirpA}\left(\frac{n}{2\mathrm{\Δ}}\right)=\exp (-j\frac{\mathrm{\π}}{4{\left(\mathrm{\Δ}\right)}^2}{n}^2\tan \frac{\mathrm{\φ}}{2}),$ | n|≤N-1, φ ≡ a pi/2, a is the conversion order;

$l_0\left(\frac{n}{\mathrm{\Δ}}\right)=\mathrm{chirpA}\left(\frac{2n}{2\mathrm{\Δ}}\right),$ $l_1\left(\frac{n}{\mathrm{\Δ}}\right)=\mathrm{chirpA}\left(\frac{2n+1}{2\mathrm{\Δ}}\right);$

3. g _iTwo heterogeneous component e of (k/ Δ) and another chirp signal chirpB (n/ (2 Δ)) _i(k/ Δ) makes convolution respectively, wherein:

$\mathrm{chirpB}\left(\frac{n}{2\mathrm{\Δ}}\right)=\exp \left(j\frac{\mathrm{\π}}{4{\left(\mathrm{\Δ}\right)}^2}{n}^2\csc \mathrm{\φ}\right),$ |n|≤2N-1；

$e_0\left(\frac{n}{\mathrm{\Δ}}\right)=\mathrm{chirpB}\left(\frac{2n}{2\mathrm{\Δ}}\right),$ $e_1\left(\frac{n}{\mathrm{\Δ}}\right)=\mathrm{chirpB}\left(\frac{2n+1}{2\mathrm{\Δ}}\right);$

4. the two-way convolution results of step in 3. intercepted the N point respectively;

The N point convolution results addition of two branch roads that 5. 4. step obtained, addition result again with the sequence l of step in 2. ₀(k/ Δ) and coefficient A _φ/ (2 Δ) multiplies each other, output end product { F ^aF} (m/ Δ), wherein | m|≤(N-1)/2, A _φ≡ exp (i π sgn (sin φ)/4+i φ/2)/| sin φ |.

Below the derivation of foregoing invention content is carried out brief description:

Suppose that input signal sequence is:

$\left|n\right|\≤\frac{N-1}{2},$ N represents that sequence length is an odd number

The normalization sampling rate is $\mathrm{\Δ}=\sqrt{N},$ Through after the interpolation, sampling rate becomes 2 Δs, and new sequence length becomes 2N-1:

|n|≤N-1

If anti-aliasing interpolation filter is h (n).

Then null value interpolation postorder is classified as:

$v\left(\frac{n}{2\mathrm{\Δ}}\right)=\left\{\begin{array}{cc}f\left(\frac{n/2}{\mathrm{\Δ}}\right),& n=0,\±2,\±4....\\ 0,& \mathrm{otherwise}\end{array}\right.$

Obtain behind the anti-aliasing filter:

$f^{\′}\left(\frac{n}{2\mathrm{\Δ}}\right)=\underset{k}{\mathrm{\Σ}}v\left(\frac{k}{2\mathrm{\Δ}}\right)h(n-k)=\underset{k}{\mathrm{\Σ}}f\left(\frac{k}{\mathrm{\Δ}}\right)h(n-2k),$

Wherein-(N-1)≤n≤N-1 ,-(N-1)/2≤k≤(N-1)/2.

If

r ₀(n)＝h(2n+1)

r ₁(n)＝h(2n)

When n=2m:

$f_1\left(\frac{m}{\mathrm{\Δ}}\right)=f^{\′}\left(\frac{2m}{2\mathrm{\Δ}}\right)=\underset{k}{\mathrm{\Σ}}f\left(\frac{k}{\mathrm{\Δ}}\right)h(2m-2k)=\underset{k}{\mathrm{\Σ}}f\left(\frac{k}{\mathrm{\Δ}}\right)r_1(m-k)$

When n=2m+1:

$f_0\left(\frac{m}{\mathrm{\Δ}}\right)=f^{\′}\left(\frac{2m+1}{2\mathrm{\Δ}}\right)=\underset{k}{\mathrm{\Σ}}f\left(\frac{k}{\mathrm{\Δ}}\right)h(2m+1-2k)=\underset{k}{\mathrm{\Σ}}f\left(\frac{k}{\mathrm{\Δ}}\right)r_0(m-k)$

Again

$\mathrm{chirpA}\left(\frac{n}{2\mathrm{\Δ}}\right)=\exp (-j\frac{\mathrm{\π}}{4{\left(\mathrm{\Δ}\right)}^2}{n}^2\tan \frac{\mathrm{\φ}}{2}),$ |n|≤N-1

$\mathrm{chirpB}\left(\frac{n}{2\mathrm{\Δ}}\right)=\exp \left(j\frac{\mathrm{\π}}{4{\left(\mathrm{\Δ}\right)}^2}{n}^2\csc \mathrm{\φ}\right),$ |n|≤2N-1

φ ≡ a pi/2 wherein, a is the conversion order.Above the heterogeneous component of two chirp sequences as follows respectively:

$e_0\left(\frac{n}{\mathrm{\Δ}}\right)=\mathrm{chirpB}\left(\frac{2n}{2\mathrm{\Δ}}\right),$ $e_1\left(\frac{n}{\mathrm{\Δ}}\right)=\mathrm{chirpB}\left(\frac{2n+1}{2\mathrm{\Δ}}\right)$

$l_0\left(\frac{n}{\mathrm{\Δ}}\right)=\mathrm{chirpA}\left(\frac{2n}{2\mathrm{\Δ}}\right),$ $l_1\left(\frac{n}{\mathrm{\Δ}}\right)=\mathrm{chirpA}\left(\frac{2n+1}{2\mathrm{\Δ}}\right)$

Then order

$g\left(\frac{n}{2\mathrm{\Δ}}\right)=\mathrm{chirpA}\left(\frac{n}{2\mathrm{\Δ}}\right)*f^{\′}\left(\frac{n}{2\mathrm{\Δ}}\right),$ |n|≤N-1

$G^a\left(\frac{m}{2\mathrm{\Δ}}\right)=\underset{k}{\mathrm{\Σ}}g\left(\frac{k}{2\mathrm{\Δ}}\right)\mathrm{chirpB}(\frac{m}{2\mathrm{\Δ}}-\frac{k}{2\mathrm{\Δ}}),$ |m|≤N-1

Have according to (5) formula

$\left\{F^{a}f\right\}\left(\frac{m}{2\mathrm{\Δ}}\right)=\frac{A_{\mathrm{\φ}}}{2\mathrm{\Δ}}*\mathrm{chirpA}\left(\frac{m}{2\mathrm{\Δ}}\right)*G^a\left(\frac{m}{2\mathrm{\Δ}}\right),$ |m|≤N-1 (6)

Each heterogeneous component substitution (6) formula abbreviation, and press dual-rate and extract { F ^aThe order that f} (m/ (2 Δ)) can obtain input signal sequence f (n/ Δ) is the N point fraction Fourier conversion result { F of a ^aF} (m/ Δ):

$\left\{F^{a}f\right\}\left(\frac{m}{\mathrm{\Δ}}\right)=\frac{A_{\mathrm{\φ}}}{2\mathrm{\Δ}}*l_0\left(\frac{m}{\mathrm{\Δ}}\right)*\left\{\underset{k^{\′}}{\mathrm{\Σ}}\right[\underset{k}{\mathrm{\Σ}}f\left(\frac{k}{\mathrm{\Δ}}\right)r_1(k^{\′}-k)]l_0\left(\frac{k^{\′}}{\mathrm{\Δ}}\right)e_0\left(\frac{m-k^{\′}}{\mathrm{\Δ}}\right)$

$+\underset{k^{\′}}{\mathrm{\Σ}}\left[\underset{k}{\mathrm{\Σ}}f\left(\frac{k}{\mathrm{\Δ}}\right)r_0(k^{\′}-k)\right]l_1\left(\frac{k^{\′}}{\mathrm{\Δ}}\right)e_1\left(\frac{m-k^{\′}}{\mathrm{\Δ}}\right)\},$ |m|≤(N-1)/2 (7)

Wherein | m|≤(N-1)/2, A _φ≡ exp (i π sgn (sin φ)/4+i φ/2)/| sin φ |.(7) formula can realize whole computation process by 5 steps introducing previously, and whole process has the characteristics of two branch road concurrent operations.

Need to prove that do not limit the type of interpolation filter in above-mentioned derivation, promptly this derivation all is suitable for any one interpolation filter.

The efficient implementation method of a kind of fraction Fourier conversion that the present invention proposes based on multisampling, its beneficial effect is:

(1) implementation method of the present invention's proposition is littler than the classical implementation method operand of primal algorithm.If interpolation filter length is 4N, the operand of two kinds of implementation methods is as shown in table 1 so, and the operand of efficient each branch road of implementation method of multisampling is all less than half of primal algorithm operand as can be seen from Table 1;

(2) implementation method of the present invention's proposition structurally has the parallel characteristics of two-way, can adopt parallel organization to improve computing velocity during realization, saves computing time, is particularly suitable for hardware realizations such as FPGA.Specifically can be referring to the FPGA implementation method example in " embodiment ";

(3) implementation method of the present invention's proposition selects for use the windowing half-band filter of low order as interpolation filter, saved half calculated amount of filtering interpolation link, and windowing can effectively prevent the oscillatory occurences of appearance when the conversion order is low, specifically sees the comparison in the accompanying drawing 2.As can be seen from the figure interpolation filter is lower at order, and under the not windowing situation, vibration has appearred in the fractional order transformation results, and error is big, and after the windowing, has then eliminated vibration.

The implementation method that table 1 the present invention proposes and the comparison of classical implementation method operand

Operand	Primal algorithm is classical to be realized	Multisampling realizes single branch road
			Multiple multiplier	12N log 2N+36N	(11Nlog 2N+22N)/2
Be added with number	24N log 2N+48N	(22Nlog 2N+23N)/2

Description of drawings

Fig. 1 is the efficient implementation method process flow diagram of multisampling;

Fig. 2 is the comparison that Hanning window half-band filter and other interpolation filter influence the fraction Fourier conversion result;

Fig. 3 is that an input signal length is the FPGA implementation structure block diagram of the inventive method of 1023;

Embodiment

Below in conjunction with accompanying drawing and FPGA embodiment summary of the invention is elaborated.

The present invention relates to the efficient implementation method of a kind of fraction Fourier conversion based on multisampling, its principle is seen formula (7), the algorithm flow chart of implementation method as shown in Figure 1, whole flow process resolves into following 5 steps finishes computing:

1. with the heterogeneous component r of input signal sequence f (n/ Δ) by interpolation filter h (n) _i(k) carry out filtering, obtain sequence f _i(m/ Δ), i=0 wherein, 1 an expression way.|n|≤(N-1)/2， $\mathrm{\Δ}=\sqrt{N}$ The normalized sampling interval of expression input signal sequence, N represents input signal sequence length, is odd number, wherein

H (n)=2w (n) sin (π n/2)/(π n), | n|≤L, 2L+1 are filter length;

w(n)＝0.5·[1+cos(nπ/L)]，|n|≤L；

r ₀(n)＝h(2n+1)，r ₁(n)＝h(2n)＝δ(n)；

2. use two heterogeneous component l of chirp signal chirpA (n/ (2 Δ)) _i(k/ Δ) respectively with the 1. filtering f as a result of respective branch of step _i(m/ Δ) multiplies each other and obtains g _i(k/ Δ), wherein

$\mathrm{chirpA}\left(\frac{n}{2\mathrm{\Δ}}\right)=\exp (-j\frac{\mathrm{\π}}{4{\left(\mathrm{\Δ}\right)}^2}{n}^2\tan \frac{\mathrm{\φ}}{2}),$ | n|≤N-1, φ ≡ a pi/2, a is the conversion order;

$l_0\left(\frac{n}{\mathrm{\Δ}}\right)=\mathrm{chirpA}\left(\frac{2n}{2\mathrm{\Δ}}\right),l_1\left(\frac{n}{\mathrm{\Δ}}\right)=\mathrm{chirpA}\left(\frac{2n+1}{2\mathrm{\Δ}}\right);$

3. g _iTwo heterogeneous component e of (k/ Δ) and another chirp signal chirpB (n/ (2 Δ)) _i(k/ Δ) makes convolution respectively, wherein

$\mathrm{chirpB}\left(\frac{n}{2\mathrm{\Δ}}\right)=\exp \left(j\frac{\mathrm{\π}}{4{\left(\mathrm{\Δ}\right)}^2}{n}^2\csc \mathrm{\φ}\right),$ |n|≤2N-1；

$e_0\left(\frac{n}{\mathrm{\Δ}}\right)=\mathrm{chirpB}\left(\frac{2n}{2\mathrm{\Δ}}\right),$ $e_1\left(\frac{n}{\mathrm{\Δ}}\right)=\mathrm{chirpB}\left(\frac{2n+1}{2\mathrm{\Δ}}\right);$

4. the two-way convolution results of step in 3. intercepted the N point respectively;

The N point convolution results addition of two branch roads that 5. 4. step obtained, addition result again with the sequence l of step in 2. ₀(k/ Δ) and coefficient A _φ/ (2 Δ) multiplies each other and exports end product { F ^aF} (m/ Δ), wherein, | m|≤(N-1)/2, A _φ≡ exp (i π sgn (sin φ)/4+i φ/2)/| sin φ |.

Provide this algorithm below in conjunction with above-mentioned 5 steps and be used for the example that FPGA realizes, the length of input signal is made as 1023 points, and interpolation filter is that the Hanning window on 123 rank partly is with the FIR wave filter.Fig. 3 is that this routine FPGA realizes theory diagram, has omitted control module among the figure, is made up of following basic module: the filtering interpolation module; The ChirpDDS module; FFT and IFFT module; Output unit; Control module.

(1) the filtering interpolation module is that the Hanning window on one 123 rank partly is with FIR wave filter h (n), owing to adopt multisampling structure, r ₁(n)=and δ (n), therefore the filtering interpolation module of second branch road can be omitted.Consider the symmetry of coefficient again, first branch road only needs the FIR wave filter on one 31 rank just can realize, wherein

$h\left(n\right)=w\left(n\right)\·\frac{2\sin \left(\frac{\mathrm{\πn}}{2}\right)}{\mathrm{\πn}},$ $w\left(n\right)=0.5\·[1+\cos \left(\frac{\mathrm{n\π}}{L}\right)],$ n＝-L，…-1，0，1，…L，L＝61

(2) the ChirpDDS module adopts traditional Chirp DDS method for designing, and being added up by phase-accumulated unit and frequency, control module is added in the unit and sine lookup table constitutes.Major function is several chirp sequences corresponding in the output formula (6): chirpA (n/2 Δ), chirpB (n/2 Δ).Two heterogeneous component l of sequence chirpA ₀(k/ Δ), l ₁(k/ Δ) is used for the 2. modulation of two branch roads of step.Two heterogeneous component e of sequence chirpB ₀(k/ Δ), e ₁(k/ Δ) is used for step convolution algorithm 3., wherein

$\mathrm{chirpA}\left(\frac{n}{2\mathrm{\Δ}}\right)=\exp (-j\frac{\mathrm{\π}}{4{\left(\mathrm{\Δ}\right)}^2}{n}^2\tan \frac{\mathrm{\φ}}{2}),$ |n|≤N-1；

$\mathrm{chirpB}\left(\frac{n}{2\mathrm{\Δ}}\right)=\exp \left(j\frac{\mathrm{\π}}{4{\left(\mathrm{\Δ}\right)}^2}{n}^2\csc \mathrm{\φ}\right),$ |n|≤2N-1；

$e_0\left(\frac{n}{\mathrm{\Δ}}\right)=\mathrm{chirpB}\left(\frac{2n}{2\mathrm{\Δ}}\right),$ $e_1\left(\frac{n}{\mathrm{\Δ}}\right)=\mathrm{chirpB}\left(\frac{2n+1}{2\mathrm{\Δ}}\right);$

$l_0\left(\frac{n}{\mathrm{\Δ}}\right)=\mathrm{chirpA}\left(\frac{2n}{2\mathrm{\Δ}}\right),$ $l_1\left(\frac{n}{\mathrm{\Δ}}\right)=\mathrm{chirpA}\left(\frac{2n+1}{2\mathrm{\Δ}}\right);$

(3) FFT module and IFFT module completing steps convolution algorithm 3., two each FFT modules of branch road are share an IFFT module at last.3 modules are carried out the length of FFT conversion by controller control, and maximum length all is 2 times of list entries length, and promptly 2N ≈ 2048;

(4) output unit is by a multiple multiplier and a buffer memory sequence l ₀The ram of (k/ Δ) and constant forms.4. and 5. corresponding step.From ram, read Chang Xulie and IFFT and export just that the convolution output multiplication obtains 1023 final fraction Fourier conversion results;

(5) control module is the state machine of a complexity, and its major function is as follows: the frequency modulation rate and the initial frequency that control the ChirpDDS output sequence according to the size of conversion order; The transform length of control FFT and IFFT module; Read-write enable signal, read/write address and other control signals of control ram.

Supposing now need be to the continuous 100 order fraction Fourier conversion that carry out of a frame N=1023 point data, whole module realizes in Xilinx xc2vp20FPGA chip, be operated in the 100M clock, the time of so whole computation process consumption is approximately 4.1ms, about this 9.5ms than the FPGA realization of classic method is fast again.

Above-described specific descriptions; purpose, technical scheme and beneficial effect to invention further describe; institute is understood that; the above only is specific embodiments of the invention; and be not intended to limit the scope of the invention; within the spirit and principles in the present invention all, any modification of being made, be equal to replacement, improvement etc., all should be included within protection scope of the present invention.

Claims (1)

1. the fraction Fourier conversion implementation method based on multisampling is characterized in that, comprises following five steps, and each step comprises the branch road of two concurrent operations:

1. with the heterogeneous component r of input signal sequence f (n/ Δ) by interpolation filter h (n) _i(k) carry out filtering, obtain sequence f _i(m/ Δ), i=0 wherein, 1 an expression way; | n|≤(N-1)/2, $\mathrm{\Δ}=\sqrt{N}$ The normalized sampling interval of expression input signal sequence, N represents input signal sequence length, is odd number, wherein:

H (n)=2w (n) sin (π n/2)/(π n), | n|≤L, 2L+1 are filter length;

w(n)＝0.5·[1+cos(nπ/L)]，|n|≤L；

r ₀(n)＝h(2n+1)，r ₁(n)＝h(2n)＝δ(n)；

2. use two heterogeneous component l of chirp signal chirpA (n/ (2 Δ)) _i(k/ Δ) respectively with the 1. filtering f as a result of respective branch of step _i(m/ Δ) multiplies each other and obtains g _i(k/ Δ), wherein:

$\mathrm{chirpA}\left(\frac{n}{2\mathrm{\Δ}}\right)=\exp (-j\frac{\mathrm{\π}}{4{\left(\mathrm{\Δ}\right)}^2}{n}^2\tan \frac{\mathrm{\φ}}{2}),$ | n|≤N-1, φ ≡ a pi/2, a is the conversion order;

$l_0\left(\frac{n}{\mathrm{\Δ}}\right)=\mathrm{chirpA}\left(\frac{2n}{2\mathrm{\Δ}}\right),$ $l_1\left(\frac{n}{\mathrm{\Δ}}\right)=\mathrm{chirpA}\left(\frac{2n+1}{2\mathrm{\Δ}}\right);$

3. g _iTwo heterogeneous component e of (k/ Δ) and another chirp signal chirpB (n/ (2 Δ)) _i(k/ Δ) makes convolution respectively, wherein:

$\mathrm{chirpB}\left(\frac{n}{2\mathrm{\Δ}}\right)=\exp \left(j\frac{\mathrm{\π}}{4{\left(\mathrm{\Δ}\right)}^2}{n}^2\csc \mathrm{\φ}\right),$ |n|≤2N-1；

$e_0\left(\frac{n}{\mathrm{\Δ}}\right)=\mathrm{chirpB}\left(\frac{2n}{2\mathrm{\Δ}}\right),$ $e_1\left(\frac{n}{\mathrm{\Δ}}\right)=\mathrm{chirpB}\left(\frac{2n+1}{2\mathrm{\Δ}}\right);$

4. the two-way convolution results of step in 3. intercepted the N point respectively;

The N point convolution results addition of two branch roads that 5. 4. step obtained, addition result again with the sequence l of step in 2. ₀(k/ Δ) and coefficient A _φ/ (2 Δ) multiplies each other, output end product { F ^aF} (m/ Δ), wherein | m|≤(N-1)/2, A _φ≡ exp (i π sgn (sin φ)/4+i φ/2)/| sin φ |.

Images