Resampling technique for arbitrary sampling rate conversion

Abstract

A resampling method based on window function for flexible sampling rate conversion in broadband frequency measurement devices is described. The resampling algorithm can satisfy the requirements of different sampling rates. The frequency responses of the filter in the resampling model based on the Farrow structure are analyzed, and the design criterion of the filter in resampling model is considered. A fractional delay filter design model based on window function method is described. A fractional delay filter matrix, which is expressed by polynomial form, is constructed. Then the expression related to subfilter coefficients is obtained and subfilter coefficients are solved for by the least square method.

Claims

1. A resampling method for flexible sampling rate conversion in a broadband frequency measurement device, the method comprising: revealing frequency responses of a filter in a resampling model; specifying design conditions of the filter; completing the resampling model, wherein: a fractional delay filter design model h.sub.p(n) is an impulse response of a fractional delay filter with a delay of p, and a group delay τ of an obtained fractional delay filter is τ+p, where τ is a group delay of an ideal-low pass filter; the fractional delay filter design model h.sub.p(n) is expressed as:
h.sub.p(n)=h.sub.d(n−p).Math.w(n−p), 0≤n≤N−1 wherein h.sub.d(n) is an impulse response of the ideal low-pass digital filter; h.sub.d(n−p) is the impulse response of the ideal low-pass digital filter shifted by a delay p, w(n) is a sequence of a window function of a length N; fractional delay filter coefficients with different group delay responses are designed by taking multiple values as p being limited as follows:
p=−0.5+bδ −0.5≤p≤0.5, wherein δ is an interval and b is an integer; a filter coefficients matrix is established by substituting p into the expressed fractional delay filter design model h.sub.p(n) to obtain multiple fractional delay filter coefficients; a fractional delay filter coefficients matrix is constructed as follows: $\begin{array}{cc}\left[\begin{array}{c}{h}_{p1}\\ h_{p2}\\ \mathrm{.Math.}\\ h_{\mathrm{pk}}\end{array}\right]=\left[\begin{array}{cccc}{h}_{p1}\left(0\right)& h_{p1}\left(1\right)& \mathrm{.Math.}& h_{p1}\left(N-1\right)\\ h_{p2}\left(0\right)& h_{p2}\left(1\right)& \mathrm{.Math.}& h_{p2}\left(N-1\right)\\ \mathrm{.Math.}& & & \\ h_{\mathrm{pk}}\left(0\right)& h_{pk}\left(1\right)& \mathrm{.Math.}& h_{pk}\left(N-1\right)\end{array}\right],& \end{array}$ wherein h.sub.pk is the filter coefficients whose delay is p.sub.k.Math.h.sub.pk=[h.sub.pk(0), h.sub.pk(1), . . . h.sub.pk(N−1)]; a polynomial form of the filter coefficients matrix is established, wherein the fractional delay filter coefficients matrix is expressed in polynomial form: $\left[\begin{array}{c}{h}_{p1}\\ h_{p2}\\ \mathrm{.Math.}\\ h_{\mathrm{pk}}\end{array}\right]=\left[\begin{array}{cccc}{p}_1& p_1^2& \mathrm{.Math.}& p_1^I\\ p_2& p_2^2& \mathrm{.Math.}& p_2^I\\ \mathrm{.Math.}& \mathrm{.Math.}& & \mathrm{.Math.}\\ p_k& p_k^2& \mathrm{.Math.}& p_k^I\end{array}\right]\left[\begin{array}{cccc}{c}_{11}& c_{12}& \mathrm{.Math.}& c_{1N}\\ c_{21}& c_{22}& \mathrm{.Math.}& c_{1N}\\ \mathrm{.Math.}& \mathrm{.Math.}& & \mathrm{.Math.}\\ c_{I1}& c_{I2}& \mathrm{.Math.}& c_{\mathrm{IN}}\end{array}\right]$ subfilter coefficients are solved by least square method in which the following
H=P⋅C is used to solve C by least square method as follows:
C=(P.sup.TP).sup.−1P.sup.TH wherein, a solution provides values for the subfilter coefficients; the values for the subfilter coefficients are used to complete the resampling model by substituting in: $y\left(k{T}_r\right)=\underset{m=-\infty }{\overset{\infty }{.Math.}}\underset{i=0}{\overset{I}{.Math.}}x\left(m{T}_s\right)c\left(\left[\frac{k{T}_r-m{T}_s}{T_r}\right],i\right){\left(\frac{m{T}_s}{T_r}-\left[\frac{m{T}_s}{T_r}\right]-0.5\right)}^i=\underset{m=-\infty }{\overset{\infty }{.Math.}}\underset{i=0}{\overset{I}{.Math.}}x\left(m{T}_s\right)c\left(\left[\frac{k{T}_r-m{T}_s}{T_r}\right],i\right){\left(p_m\right)}^i,$ wherein x(mT.sub.s) is input samples at sampling cycles T.sub.s, and T.sub.r is a resampling cycle, c(n,i) is the subfilter coefficient c.sub.n,i, and [*] indicates the largest integer not exceeding (*); wherein p.sub.m is associated with m, and $y\left(k{T}_r\right)=\underset{m=-\infty }{\overset{\infty }{.Math.}}\underset{i=0}{\overset{I}{.Math.}}{\hat{x}}_i\left(m{T}_s\right)c\left(\left[\frac{k{T}_r-m{T}_s}{T_r}\right],i\right)$ is implemented by the Transposed Farrow structure; where {circumflex over (x)}.sub.i(mT.sub.s) is a representation of a sum of x(mT.sub.s) and p.sub.m.sup.i; and using the completed resampling model in the broadband frequency measurement device to achieve a predetermined sampling rate conversion and provide conditions for measuring different component signals.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) In order to facilitate a fuller understanding of the present invention, reference is now made to the attached drawings. The drawings should not be construed as limiting the present invention, but are intended only to illustrate different aspects and embodiments of the invention.

(2) FIG. 1 shows a flowchart for arbitrary sampling rate conversion.

(3) FIG. 2 shows exemplary results and spectrum of sampling from analog filters.

(4) FIG. 3 shows the flowchart of an exemplary resampling algorithm based on window function.

(5) FIG. 4 shows the magnitude error of an exemplary VFD filter.

(6) FIG. 5 shows the group delay error of an exemplary VFD filter.

(7) FIG. 6 shows the magnitude responses of the corresponding exemplary VFD filter.

(8) FIG. 7 shows the magnitude error of the designed exemplary VFD filter.

(9) FIG. 8 shows the group delay error of the designed exemplary VFD filter.

(10) FIG. 9 shows the magnitude responses of the exemplary VFD filter.

(11) FIG. 10 illustrates exemplary hardware components for a server.

DETAILED DESCRIPTION

(12) Exemplary embodiments of the invention will now be described in order to illustrate various features of the invention. The embodiments described herein are not intended to be limiting as to the scope of the invention, but rather are intended to provide examples of the components, use, and operation of the invention.

I. Introduction

(13) With the increasing development of renewables and DC transmission, a large number of electronic devices are applied to power systems, which complicate the calculation of the electrical quantity. Therefore, the measurement devices are required to monitor the components within a wide frequency bandwidth, instead of just the fundamental one, such as the broadband frequency measurement devices in China and the Universal Grid Analyzer in the United States. Taking the broadband frequency measurement device as an example, besides the fundamental component, it needs to measure harmonics up to 50.sup.th order and inter-harmonics from 0 to 300 Hz.

(14) To measure the signals distributed within such a wide range, it is necessary to use different sampling frequencies for the calculation of different frequency components to decrease the computational burden. Therefore, the resampling process after the original sampling is critical. The resampling algorithm applied in the measurement devices requires the function of arbitrary sampling rate conversion and the prevention of frequency aliasing. In recent decades, certain resampling technologies have been widely used in radio, audio, digital signal processing, asynchronous sampling, and frequency estimation, among others.

(15) The existing methods for the implementation of resampling algorithms can be categorized into three types. The first is a direct method of decimation, interpolation and combination, the second is based on a polyphase filter structure, and the third is based on polynomial interpolation. The first method is simple. However, once the sampling rate is changed, the filter needs to be redesigned; therefore, this method has poor flexibility. The second method uses a polyphase filter-bank to achieve synchronization and interpolation. The polyphase filter structure needs to design multiple different filter units to achieve the parallel filtering function. Compared with the first method, it can realize arbitrary rational sampling rate conversion; however, it requires multiple filter units, which creates a large calculational burden.

(16) The principle of the third type of method is to implement interpolation via polynomial function, which can be further divided into Lagrange polynomials, spline interpolation polynomials and methods based on Farrow structure. The filters designed based on the Lagrange and spline interpolation polynomials have stable magnitude responses in the low-frequency band, but poor magnitude responses in the high-frequency band; thus, the performance of these filters is not great. An interpolation filter has been proposed, the impulse response of which is composed of a piecewise polynomial form that has good flexibility. When resampling by this method, the Farrow structure subfilter coefficients do not need to be changed, and the interpolation filter coefficients can be calculated according to different fractional delay parameters. This method can achieve efficient calculation and arbitrary sampling rate conversion, and has low calculational burden.

(17) The key to design the Farrow structure interpolation filter is the determination of the subfilter coefficients. The existing coefficients design methods include weighted least squares (WLS), minimax methods, and intelligent algorithms. Among them, intelligent algorithms represented by particle swarm optimization (PSO) are only applicable to the cases of low-dimension coefficients; when the dimension of coefficients is high, the algorithm cannot be further optimized. The minimax method and WLS method are both designed to minimize the design error. During the design process, the matrix inversion and integration operations increase the calculation complexity, resulting in a long design time and the complicated adjustment of weight parameters. In addition, a matrix-based conjugate gradient (CG) algorithm has been considered, and this method can effectively solve the problems of large amounts of calculation and weight adjustment in optimal algorithms. However, the focus of these methods is on the design of the filter in the passband range, and they do not consider the transition band or stopband. Therefore, whether the performance of the filters designed by these methods can meet the actual requirements needs to be further investigated.

(18) In this disclosure, a resampling algorithm based on window function is described, which can realize arbitrary sampling rate conversion and adjust the filter characteristics according to the demand. According to the resampling model based on Farrow structure, the change principles of the filter characteristics are investigated, namely the larger the filter order, the shorter the frequency band in the spectrum, and the design criterion that requires the stopband cutoff frequency of the filter in the resampling model less than π is described. Then a method based on window function is described to obtain multiple fractional delay filter coefficients, subfilter coefficients matrix is solved by least square method, and the characteristics of window spectrum are described. Simulation results show that resampling based on the presented method can achieve arbitrary sampling rate conversion, and resampling satisfies the Nyquist theorem.

II. Characteristic Analysis of Resampling Based on Farrow Structure

(19) A. Principle of Resampling Based on Farrow Structure

(20) The process of arbitrary sampling rate conversion can be viewed as first reconstructing the original signal with the sampling points, then resampling the reconstructed signal at the required time and sampling rate. FIG. 1 shows a flowchart for arbitrary sampling rate conversion.

(21) The expression of resampling sample y(kT.sub.r) is

(22) $\begin{array}{cc}y\left(k{T}_r\right)=\underset{m=-\infty }{\overset{\infty }{.Math.}}x\left(m{T}_s\right)h_a\left(k{T}_r-m{T}_s\right),& \left(1\right)\end{array}$
where k is the number of resampling sequence, T.sub.r is the resampling cycle, m is the number of sampling sequence, T.sub.s is the sampling cycle, x(mT.sub.s) is the sampling sequence, h.sub.a(t) is the impulse response of the ideal analog low-pass filter, and h.sub.a(kT.sub.r−mT.sub.s) represents the sample of h.sub.a(t). Hence, the resampling result (1) can be considered the convolution of x(mT.sub.s) with the impulse response h.sub.d(m) of digital filter.

(23) To ensure that frequency aliasing does not occur during resampling, the frequency responses H.sub.d(ω) of the digital filter h.sub.d(m) should be

(24) $\begin{array}{cc}{H}_d\left(\omega \right)=\{\begin{array}{cc}1,& .Math.\omega .Math.<f_r/2,\\ 0,& .Math.\omega .Math.\ge f_r/2,\end{array}& \left(2\right)\end{array}$
where ω is the normalized angular frequency and f.sub.r is the resampling frequency. However, it is impossible to realize the impulse response h.sub.a(t) of the ideal analog low-pass filter in practical applications.

(25) Simple piecewise polynomials are generally used instead of the impulse response h.sub.a(t). Given the degree I of the polynomial and the length Δ, h.sub.a(t) can be divided into equal-length piecewise polynomials

(26) $\begin{array}{cc}\begin{array}{cc}{h}_a\left(t\right)=\underset{n=0}{\overset{N-1}{.Math.}}{a}_n\left(t-n\Delta \right),& 0\le t\le N\Delta ,\end{array}& \left(3\right)\end{array}$
where N is the number of polynomials and a.sub.n(t) is the n-th segment polynomial. a.sub.n(t) is expressed as

(27) $\begin{array}{cc}{a}_n\left(t\right)=\{\begin{array}{cc}\underset{i=0}{\overset{I}{.Math.}}c\left(n,i\right){\left(\frac{t}{\Delta }-0\mathrm{.5}\right)}^i,& 0\le t\le \Delta ,\\ 0,& \mathrm{esle}t,\end{array}& \left(4\right)\end{array}$
where c(n,i) is the i-th coefficient of the n-th segment polynomial, and t/Δ−0.5 is the polynomial variable with a range of [−0.5, 0.5]. After taking n=[t/Δ], and substituting (4) into (3),

(28) $\begin{array}{cc}\begin{array}{cc}{h}_a\left(t\right)=\underset{i=0}{\overset{I}{.Math.}}c\left(\left[\frac{t}{\Delta }\right],i\right){\left(\frac{t}{\Delta }-\left[\frac{t}{\Delta }\right]-0.5\right)}^i,& 0\le t\le N\Delta ,\end{array}& \left(5\right)\end{array}$
where [⋅] indicates the largest integer not exceeding (⋅). Equation (5) is the expression of replacing h.sub.a(t) with piecewise polynomials.

(29) After substituting (5) into (1) and taking Δ=T.sub.s, the expression becomes

(30) $\begin{array}{cc}\begin{array}{c}y\left(k{T}_r\right)=\underset{m=-\infty }{\overset{\infty }{.Math.}}\underset{i=0}{\overset{I}{.Math.}}x\left(m{T}_s\right)c\left(\left[\frac{k{T}_r-m{T}_s}{T_s}\right],i\right){\left(\frac{k{T}_r}{T_s}-\left[\frac{k{T}_r}{T_s}\right]-0\mathrm{.5}\right)}^i\\ =\underset{m=-\infty }{\overset{\infty }{.Math.}}\underset{i=0}{\overset{I}{.Math.}}x\left(m{T}_s\right)c\left(\left[\frac{k{T}_r-m{T}_s}{T_s}\right],i\right){\left(p_k\right)}^i,\end{array}& \left(6\right)\end{array}$
where p.sub.k is the delay between the input sample and the output sample. As p.sub.k has nothing to do with m, (6) can be simplified to

(31) $\begin{array}{cc}y\left({\mathrm{kT}}_r\right)=\underset{i=0}{\overset{I}{.Math.}}{\tilde{x}}_i\left({\mathrm{mT}}_r\right){\left(p_k\right)}^i,& \left(7\right)\end{array}$
where {tilde over (x)}.sub.i(mT.sub.r) is the representation of the sum of x(mT.sub.s) and c([(kT.sub.r−mT.sub.s)/T.sub.s,i]). Equation (7) can be implemented by the Farrow structure.

(32) After substituting (5) into (1) and taking Δ=T.sub.r, the expression becomes

(33) $\begin{array}{cc}\begin{array}{c}y\left(k{T}_r\right)=\underset{m=-\infty }{\overset{\infty }{.Math.}}\underset{i=0}{\overset{I}{.Math.}}x\left(m{T}_s\right)c\left(\left[\frac{k{T}_r-m{T}_s}{T_r}\right],i\right){\left(\frac{-{\mathrm{mT}}_s}{T_r}-\left[-\frac{{\mathrm{mT}}_s}{T_r}\right]-0\mathrm{.5}\right)}^i\\ =\underset{m=-\infty }{\overset{\infty }{.Math.}}\underset{i=0}{\overset{I}{.Math.}}x\left(m{T}_s\right)c\left(\left[\frac{k{T}_r-m{T}_s}{T_r}\right],i\right){\left(p_m\right)}^i.\end{array}& \left(8\right)\end{array}$

(34) As p.sub.m is associated with m, (8) can be simplified to

(35) $\begin{array}{cc}y\left(k{T}_r\right)=\underset{m=-\infty }{\overset{\infty }{.Math.}}\underset{i=0}{\overset{I}{.Math.}}{\hat{x}}_i\left(m{T}_s\right)c\left(\left[\frac{k{T}_r-m{T}_s}{T_r}\right],i\right),& \left(9\right)\end{array}$
where {tilde over (x)}.sub.i(mT.sub.s) is the representation of the sum of x(mT.sub.s) and p.sub.m.sup.i. Equation (9) can be implemented by the Transposed Farrow structure. Equations (6) and (8) can both be denoted as the resampling model based on Farrow structure. Because the segment length is different, the filter in (6) has better de-mirror properties, and (6) is suitable for interpolation. The filter in (8) has relatively better de-aliasing properties, and (8) is suitable for decimation. In Farrow structure, c(n,i) are called the subfilter coefficients.

(36) Because there are delay variables p.sub.k and p.sub.m in (6) and (8), the group delay responses of the filter in the resampling model can be changed with adjustable delay.

(37) In resampling models (6) and (8), the key is to determine the subfilter coefficients c(n,i). Additionally, to ensure that resampling meets the Nyquist theorem, the digital filtering fitted by c(n,i) is required to filter out signals above f.sub.r/2.

(38) B. Frequency Responses Analysis of the Filter in Resampling Model

(39) Because the resampling model must have the function of arbitrary sampling rate conversion, the frequency responses of the filter in the model should change with the resampling frequency f.sub.r as shown in (2).

(40) Because (3) specifies the range of t in h.sub.d(t), the range of m in (1) is specified correspondingly as
0≤kT.sub.r−mT.sub.s≤NΔ,  (10)
and the following can be obtained:

(41) 0$\begin{array}{cc}\frac{{\mathrm{kT}}_r-N\Delta }{T_s}\le m\le \frac{{\mathrm{kT}}_r}{T_s}.& \left(11\right)\end{array}$

(42) As m is an integer, the order K of the digital filter is determined as
K=[NΔ/T.sub.s]−1.  (12)

(43) When Δ=T.sub.s, K is a fixed value. Equation (6) is suitable for interpolation, that is, the case in which T.sub.s≥T.sub.r. If the digital filter satisfies (2) when T.sub.s=T.sub.r, then resampling can satisfy the Nyquist theorem when T.sub.s>T.sub.r.

(44) When Δ=T.sub.r, K varies with T.sub.r/T.sub.s. Equation (8) is suitable for decimation, that is, the case in which T.sub.s≤T.sub.r. If T.sub.s is specified as a fixed value, then the frequency responses of the filter need to change with T.sub.r to satisfy (2). Therefore, the frequency responses of the filter in (8) are analyzed as follows.

(45) According to (5), h.sub.p(n) is defined as the impulse response of the digital filter corresponding to h.sub.a(t). h.sub.p(n) can be simply expressed as

(46) $\begin{array}{cc}\begin{array}{cc}{h}_p\left(n\right)=\underset{i=0}{\overset{I}{.Math.}}c\left(n,i\right)p^i,& -0.5\le p\le 0.5.\end{array}& \left(13\right)\end{array}$

(47) Because the group delay responses of h.sub.p(n) vary with delay p, (13) represents the impulse response of variable fractional delay (VFD) filter. When p is fixed, h.sub.p(n) is the impulse response of a fractional delay (FD) filter.

(48) Equation (1) shows that h.sub.p(n) can be regarded as the result of sampling h.sub.a(t). Take the sampling period as T.sub.1, then
h.sub.p(n)=h.sub.a(t)|.sub.t=nT.sub.1.  (14)

(49) To determine the relationship between h.sub.p(n) and h.sub.a(t) in the frequency domain, an ideal sample signal h.sub.a′(t) is introduced. The expression of h.sub.a′(t) is

(50) $\begin{array}{cc}{h}_a^{\prime }\left(t\right)=\underset{n=-\infty }{\overset{\infty }{.Math.}}{h}_a\left(t\right)\delta \left(t-n{T}_1\right),& \left(15\right)\end{array}$
where δ(t) is the impulse function. Take the Fourier transforms corresponding to h.sub.a(t), h.sub.p(n), and h.sub.a′(t) as H.sub.a(Ω), H.sub.p(ω) and H.sub.a′(Ω), there are

(51) $\begin{array}{cc}{H}_p\left(\omega \right)=\underset{n=-\infty }{\overset{\infty }{.Math.}}h\left(n\right)e^{-j\omega n},& \left(16\right)\end{array}$
and

(52) $\begin{array}{cc}\begin{array}{c}{H}_a^{\prime }\left(\Omega \right)={\int }_{-\infty }^{\infty }{h}_a^{\prime }\left(t\right)e^{-j\Omega t}dt={\int }_{-\infty }^{\infty }\underset{n=-\infty }{\overset{\infty }{.Math.}}{h}_a^{\prime }\left(t\right)\delta \left(t-n{T}_1\right)e^{-j\Omega t}\mathrm{dt}\\ =\underset{n=-\infty }{\overset{\infty }{.Math.}}{h}_a\left(n{T}_1\right)e^{-j\Omega T_{1}n}=\underset{n=-\infty }{\overset{\infty }{.Math.}}{h}_a\left(n{T}_1\right)e^{-j\omega n},\end{array}& \left(17\right)\end{array}$
where Ω is the analog angular frequency, Ω=ωT.sub.1. Therefore, the relationship between H.sub.p(ω) and H.sub.a′(Ω) is
H.sub.p(ω)=H.sub.a′(Ω).  (18)

(53) The relationship between H.sub.a′(Ω) and H.sub.a(Ω) is

(54) $\begin{array}{cc}{H}_a^{\prime }\left(\Omega \right)=\frac{1}{2\pi }\left[\mathrm{FT}\left[h_a\left(t\right)\right]*FT\left[\delta \left(t-n{T}_1\right)\right]\right]=\frac{1}{T_1}\underset{n=-\infty }{\overset{\infty }{.Math.}}{H}_a\left(\Omega -\frac{2\pi n}{T_1}\right),& \left(19\right)\end{array}$
where FT[⋅] is Fourier transform. By combining (18) and (19), the relationship between H.sub.p(ω) and H.sub.a(Ω) is obtained:

(55) $\begin{array}{cc}{H}_p\left(\omega \right)=H_a^{\prime }\left(\Omega \right)=\frac{1}{T_1}\underset{n=-\infty }{\overset{\infty }{.Math.}}{H}_a\left(\Omega -\frac{2\pi n}{T_1}\right)=\frac{1}{T_1}\underset{n=-\infty }{\overset{\infty }{.Math.}}{H}_a\left(\frac{\omega -2\pi n}{T_1}\right).& \left(20\right)\end{array}$

(56) Equation (20) shows that H.sub.p(ω) is obtained by expanding H.sub.a(Ω) by T.sub.1 times, multiplying it by 1/T.sub.1, and then shifting it by 2πn/T.sub.1.

(57) To reflect the frequency domain characteristics of filters with different orders, another sampling period T.sub.2 (T.sub.2>T.sub.1) is used to sample h.sub.a(t) to obtain h.sub.p′(n′). Its Fourier transform form is H.sub.p′(ω′). According to (20), the expression of H.sub.p′(ω′) is

(58) $\begin{array}{cc}{H}_p^{\prime }\left({\omega }^{\prime }\right)=\frac{1}{T_2}\underset{n^{\prime }=-\infty }{\overset{\infty }{.Math.}}{H}_a\left(\Omega -\frac{2\pi n^{\prime }}{T_2}\right)=\frac{1}{T_2}\underset{n^{\prime }=-\infty }{\overset{\infty }{.Math.}}{H}_a\left(\frac{{\omega }^{\prime }-2\pi n^{\prime }}{T_2}\right).& \left(21\right)\end{array}$
where ω′=ΩT.sub.2. By taking D=T.sub.2/T.sub.1, then ω′=Dω. By contrasting (20) and (21), it can be found that the larger the sampling period, the wider the frequency band.

(59) FIG. 2 shows exemplary results and spectrum of sampling from analog filters. In this example, diagrams of h.sub.a(t), h.sub.p(n), h.sub.p′(n′) and H.sub.a(Ω), H.sub.p(ω), H.sub.p′(ω′) are displayed. FIG. 2 indicates that different resampling period T.sub.r further affect the frequency responses of the filter by affecting the filter order; the larger the filter order, the shorter the frequency band in the spectrum.

(60) According to (20) and FIG. 2, H.sub.p(ω) is obtained by expanding the spectrum of H.sub.a(Ω) with a period of 2π. After specifying that T.sub.s is fixed, when T.sub.s=T.sub.r, if the stopband cutoff frequency (ω.sub.s) of the digital filter is controlled at ω.sub.s<π, resampling satisfy the Nyquist theorem. Then for T.sub.r=DT.sub.s (D>1), the Nyquist theorem requires ω.sub.s<π/D, and the cutoff frequency of the digital filter will be less than π/D according to (21). Therefore, resampling can also satisfy the sampling theorem when T.sub.r=DT.sub.s.

(61) Thus, the design criterion requires that in resampling models (6) and (8), the stopband cutoff frequency ω.sub.s of the filter should be ω.sub.s<π when T.sub.s=T.sub.r, then for cases in which T.sub.r≠T.sub.s, frequency aliasing does not occur in resampling. Therefore, for the determination of the subfilter coefficients c(n,i) in (6) and (8), only the situation in which T.sub.s=T.sub.r needs to be considered.

III. Resampling Based on Window Function Design

(62) The key to the resampling model is the solution of the subfilter coefficients. In this section, a method based on window function is described to design subfilter coefficients. To ensure the characteristics of the filter, the window function is then discussed.

(63) A. The Solution of Subfilter Coefficients in Farrow Structure

(64) The analysis in the previous section demonstrated that the design of the subfilter coefficients only needs to consider the case in which T.sub.s=T.sub.r in the resampling model; therefore, the filter order is specified as N. The frequency response corresponding to h.sub.p(n) is

(65) $\begin{array}{cc}{H}_p\left(\omega ,p\right)=\underset{n=0}{\overset{N-1}{.Math.}}{h}_p\left(n\right)e^{-j\omega n}=\underset{n=0}{\overset{N-1}{.Math.}}\underset{i=0}{\overset{I}{.Math.}}c\left(n,i\right)p^i{e}^{-j\omega n}.& \left(22\right)\end{array}$

(66) According to the requirement of (2), the desired variable frequency response is
H.sub.d(ω,p)=e.sup.−jω((N-1)/2+p),ω∈[0,απ],  (23)
where 0<α<1, [0, απ] is the passband of the filter. The traditional methods for the design of the subfilter coefficients c(n,i) are taking c(n,i) as variable, and taking the difference between H.sub.p(ω, p) and H.sub.d(ω, p) as the objective function. The objective function is then optimized, and the coefficients are finally solved. Hence, the traditional methods only focus on the difference of the passband without considering the stopband, and there are defects in these methods.

(67) Unlike traditional methods, subfilter coefficients are designed from the perspective of designing fractional delay filter in advance. In this disclosure, the window function method is used to design a fractional delay filter. The expression of designing the impulse response of the linear-phase FIR filter via window function method is
h.sub.l(n)=h.sub.d(n).Math.w(n), 0≤n≤N−1,  (24)
where w(n) is a sequence of the window function of length N, Kaiser window is selected in this disclosure, and h.sub.d(n) is the impulse response of the ideal low-pass digital filter. The expression of h.sub.d(n) is

(68) $\begin{array}{cc}{h}_d\left(n\right)=\{\begin{array}{cc}\frac{\sin \left[{\omega }_c\left(n-\tau \right)\right]}{\pi \left(n-\tau \right)},& n\ne \tau ,\\ \frac{{\omega }_c}{\pi },& n=\tau ,\end{array}& \left(25\right)\end{array}$
where τ and ω.sub.c are the group delay and the cutoff frequency of the ideal low-pass digital filter. The relationship between ω.sub.c, the passband cutoff frequency ω.sub.p and the stopband cutoff frequency ω.sub.s of the designed filter is
ω.sub.c=(ω.sub.p+ω.sub.s)/2.  (26)

(69) Thus, ω.sub.c can indirectly determine ω.sub.s. The group delay τ of the filter obtained by (24) which can be considered a fractional delay filter without delay, is (N−1)/2. To determine the impulse response h.sub.p(n) of the fractional delay filter, the sequence n is shifted by −p. Equation (24) then becomes
h.sub.p(n)=h.sub.d(n−p).Math.w(n−p), 0≤n≤N−1.  (27)

(70) Equation (27) can be seen as the group delay τ of h.sub.l(n) shifted by p. Therefore, the group delay of the obtained fractional delay filter is τ+p, and h.sub.p(n) is the impulse response of a fractional delay filter with a delay of p. The fractional delay filter obtained by window function method can be regarded as a result in which the magnitude responses of the FIR linear phase filter are unchanged, but the group delay responses are shifted.

(71) Because the fractional delay filter can be designed in advance, it is much easier to design the subfilter coefficients c(n,i). In the interval [−0.5,0.5], take multiple values as p as follows:
p=−0.5+bδ −0.5≤p≤0.5,  (28)
where δ is an interval and b is an integer. Substitute p into (27) to obtain multiple fractional delay filter coefficients. A fractional delay filter coefficients matrix can then be constructed:

(72) 0$\begin{array}{cc}\left[\begin{array}{c}{h}_{p1}\\ h_{p2}\\ \mathrm{.Math.}\\ h_{\mathrm{pk}}\end{array}\right]=\left[\begin{array}{cccc}{h}_{p1}\left(0\right)& h_{p1}\left(1\right)& \mathrm{.Math.}& h_{p1}\left(N-1\right)\\ h_{p2}\left(0\right)& h_{p2}\left(1\right)& \mathrm{.Math.}& h_{p2}\left(N-1\right)\\ \mathrm{.Math.}& & & \\ h_{\mathrm{pk}}\left(0\right)& h_{\mathrm{pk}}\left(1\right)& \mathrm{.Math.}& h_{\mathrm{pk}}\left(N-1\right)\end{array}\right],& \left(29\right)\end{array}$
where h.sub.p.sub.k is the filter coefficients whose delay is p.sub.k. h.sub.p.sub.k=[h.sub.p.sub.k(0), h.sub.p.sub.k(1), . . . , h.sub.p.sub.k(N−1)]. By combining (29) and (13), the coefficients matrix can be expressed in polynomial form:

(73) $\begin{array}{cc}\left[\begin{array}{c}{h}_{p1}\\ h_{p2}\\ \mathrm{.Math.}\\ h_{pk}\end{array}\right]=\left[\begin{array}{cccc}{p}_1& p_1^2& \mathrm{.Math.}& p_1^I\\ p_2& p_2^2& \mathrm{.Math.}& p_2^I\\ \mathrm{.Math.}& \mathrm{.Math.}& & \mathrm{.Math.}\\ p_k& p_k^2& \mathrm{.Math.}& p_k^I\end{array}\right]\left[\begin{array}{cccc}{c}_{11}& c_{12}& \mathrm{.Math.}& c_{1N}\\ c_{21}& c_{22}& \mathrm{.Math.}& c_{1N}\\ \mathrm{.Math.}& \mathrm{.Math.}& & \mathrm{.Math.}\\ c_{I1}& c_{I2}& \mathrm{.Math.}& c_{\mathrm{IN}}\end{array}\right].& \left(30\right)\end{array}$

(74) Then (30) can be simply to be
H=P□C.  (31)

(75) Equation (31) can then be used to solve C by least square method as follows:
C=(P.sup.TP).sup.−1P.sup.TH.  (32)

(76) By substituting the transpose of the solution result C.sup.T into (6) or (8), the entire model can be completely established.

(77) B. The Window Function Selection

(78) Equation (27) shows that, during the process of designing h.sub.p(n), the window function w(n) can be adjusted to change h.sub.p(n). The window function can determine the difference between (22) and (23), and adjust the stopband attenuation of the designed filter. Therefore, the anti-aliasing effect of the designed filter is affected by the window function.

(79) Different window functions have different truncation effects. Taking the Hanning and Blackman window as an example, there will be a slightly large leakage factor when truncating; thus, it is difficult to ensure the passband stability of the designed filter. As Kaiser window can adjust its frequency responses by modifying its parameters, Kaiser window with strong adaptability for design is selected in this disclosure.

(80) The expression of the Kaiser window in the time domain is

(81) $\begin{array}{cc}w\left(n\right)=\frac{I_0\left(\beta \sqrt{1-{\left(1-\frac{2n}{N-1}\right)}^2}\right)}{I_0\left(\beta \right)},& \left(33\right)\end{array}$
where I.sub.0(⋅) is the zeroth-order modified Bessel function and β is an optional parameter. I.sub.0(⋅) can be computed as

(82) $\begin{array}{cc}{I}_0\left(x\right)=\underset{k=0}{\overset{\infty }{.Math.}}{\left[\frac{1}{k!}{\left(\frac{x}{2}\right)}^k\right]}^2,& \left(34\right)\end{array}$

(83) This infinite series is usually approximated by the sum of the first 25 terms. In (33), the characteristics of Kaiser window are different for different β, and β can be determined by the desired stopband attenuation A.sub.s.

(84) $\begin{array}{cc}\beta =\{\begin{array}{cc}0.1102\left(A_s-8.7\right),& A_s\ge 50\mathrm{dB}\\ 0.5584{\left(A_s-2.1\right)}^{0.4}+0.07886\left(A_s-21\right),& 50\mathrm{dB}>A_s>21\mathrm{dB}\\ 0,& A_s\le 21\mathrm{dB}\end{array}.& \left(35\right)\end{array}$

(85) Given A.sub.s and the order of Kaiser window N, the transition bandwidth of the filter Aw can be calculated as

(86) $\begin{array}{cc}\Delta \omega =\frac{7.95-A_s}{2.285\left(N-1\right)},& \left(36\right)\end{array}$

(87) Taking N=67 as an example, the leakage factor, relative sidelobe attenuation and mainlobe width of Kaiser window spectrum for different β are calculated, and the results are exhibited in Table I.

(88) TABLE-US-00001 TABLE I KAISER WINDOW CHARACTERISTICS FOR DIFFERENT PARAMETERS Leakage Factor Relative sidelobe Mainlobe β (%) attenuation (dB) width (π) 1 6.61  −14.7 0.025 5 0.02  −37.2 0.039 10 0.00  −74.1 0.053 15 0.00 −113.8 0.063 20 0.00 −155.1 0.072

(89) It can be found that when β is small, there will be leakage factors in the window spectrum. The larger the β, the larger the sidelobe attenuation and mainlobe width of the window spectrum. When β is sufficiently large, the passband of the designed filter becomes narrower. Thus, Kaiser window can design different frequency responses according to the required stopband attenuation and transition bandwidth.

(90) C. Exemplary Flowchart

(91) FIG. 3 shows the flowchart of an exemplary resampling algorithm based on window function. In this example, the presented method is simple in operation and the matrix operation is only reflected in least square method.

(92) In step 310, a fractional delay filter design model based on window function is established. In particular, the expression of designing the impulse response of the linear-phase FIR filter via window function method is
h.sub.l(n)=h.sub.d(n).Math.w(n), 0≤n≤N−1,  (37)
where w(n) is a sequence of the window function of length N, Kaiser window is selected in this disclosure, and h.sub.d(n) is the impulse response of the ideal low-pass digital filter. The expression of h.sub.d(n) is

(93) $\begin{array}{cc}{h}_d\left(n\right)=\{\begin{array}{cc}\frac{\sin \left[{\omega }_c\left(n-\tau \right)\right]}{\pi \left(n-\tau \right)},& n\ne \tau ,\\ \frac{{\omega }_c}{\pi },& n=\tau ,\end{array}& \left(38\right)\end{array}$
where τ and ω.sub.c are the group delay and the cutoff frequency of the ideal low-pass digital filter. The relationship between ω.sub.c, the passband cutoff frequency ω.sub.p and the stopband cutoff frequency ω.sub.s of the designed filter is
ω.sub.c=(ω.sub.p+ω.sub.s)/2.  (39)

(94) Thus, ω.sub.c can indirectly determine ω.sub.s. The group delay τ of the filter obtained by (24) which can be considered a fractional delay filter without delay, is (N−1)/2. To determine the impulse response h.sub.p(n) of the fractional delay filter, the sequence n is shifted by −p. Equation (24) then becomes
h.sub.p(n)=h.sub.d(n−p).Math.w(n−p), 0≤n≤N−1.  (40)
Equation (27) can be seen as the group delay τ of h.sub.l(n) shifted by p. Therefore, the group delay of the obtained fractional delay filter is τ+p, and h.sub.p(n) is the impulse response of a fractional delay filter with a delay of p.

(95) In step 320, according to the model, fractional delay filter coefficients with different group delay responses are designed, and a filter coefficients matrix is established. In particular, because the fractional delay filter can be designed in advance, it is much easier to design the subfilter coefficients c(n,i). In the interval [−0.5,0.5], take multiple values asp as follows:
p=−0.5+bδ−0.5≤p≤0.5,  (41)
where δ is an interval and b is an integer. Substitute p into (27) to obtain multiple fractional delay filter coefficients. A fractional delay filter coefficients matrix can then be constructed:

(96) $\begin{array}{cc}\left[\begin{array}{c}{h}_{p1}\\ h_{p2}\\ \mathrm{.Math.}\\ h_{\mathrm{pk}}\end{array}\right]=\left[\begin{array}{cccc}{h}_{p1}\left(0\right)& h_{p1}\left(1\right)& \mathrm{.Math.}& h_{p1}\left(N-1\right)\\ h_{p2}\left(0\right)& h_{p2}\left(1\right)& \mathrm{.Math.}& h_{p2}\left(N-1\right)\\ \mathrm{.Math.}& & & \\ h_{\mathrm{pk}}\left(0\right)& h_{\mathrm{pk}}\left(1\right)& \mathrm{.Math.}& h_{\mathrm{pk}}\left(N-1\right)\end{array}\right],& \left(42\right)\end{array}$
where h.sub.p.sub.k is the filter coefficients whose delay is p.sub.k. h.sub.p.sub.k=[h.sub.p.sub.k(0), h.sub.p.sub.k(1), . . . , h.sub.p.sub.k(N−1)].

(97) Next, in step 330, the polynomial form of the filter coefficients matrix is established, and subfilter coefficients are solved by least square method. In particular, by combining (29) and (13), the coefficients matrix can be expressed in polynomial form:

(98) $\begin{array}{cc}\left[\begin{array}{c}{h}_{p1}\\ h_{p2}\\ \mathrm{.Math.}\\ h_{pk}\end{array}\right]=\left[\begin{array}{cccc}{p}_1& p_1^2& \mathrm{.Math.}& p_1^I\\ p_2& p_2^2& \mathrm{.Math.}& p_2^I\\ \mathrm{.Math.}& \mathrm{.Math.}& & \mathrm{.Math.}\\ p_k& p_k^2& \mathrm{.Math.}& p_k^I\end{array}\right]\left[\begin{array}{cccc}{c}_{11}& c_{12}& \mathrm{.Math.}& c_{1N}\\ c_{21}& c_{22}& \mathrm{.Math.}& c_{1N}\\ \mathrm{.Math.}& \mathrm{.Math.}& & \mathrm{.Math.}\\ c_{I1}& c_{I2}& \mathrm{.Math.}& c_{\mathrm{IN}}\end{array}\right].& \left(43\right)\end{array}$

(99) Then (30) can be simply to be
H=P□C.  (44)

(100) Equation (31) can then be used to solve C by least square method as follows:
C=(P.sup.TP).sup.−1P.sup.TH.  (45)
The solution will provide the values for the coefficient matrix.

(101) Next, in step 340, the coefficients solved in step 330 can be used to complete the resampling model. In particular, the coefficients can be substituted in equations (6) or (8). Equation (6) is:

(102) $\begin{array}{cc}\begin{array}{c}y\left(k{T}_r\right)=\underset{m=-\infty }{\overset{\infty }{.Math.}}\underset{i=0}{\overset{I}{.Math.}}x\left(m{T}_s\right)c\left(\left[\frac{k{T}_r-m{T}_s}{T_s}\right],i\right){\left(\frac{k{T}_r}{T_s}-\left[\frac{k{T}_r}{T_s}\right]-0\mathrm{.5}\right)}^i\\ =\underset{m=-\infty }{\overset{\infty }{.Math.}}\underset{i=0}{\overset{I}{.Math.}}x\left(m{T}_s\right)c\left(\left[\frac{k{T}_r-m{T}_s}{T_s}\right],i\right){\left(p_k\right)}^i,\end{array}& \left(46\right)\end{array}$
where p.sub.k is the delay between the input sample and the output sample. As p.sub.k has nothing to do with m, (6) can be simplified to

(103) 0$\begin{array}{cc}y\left({\mathrm{kT}}_r\right)=\underset{i=0}{\overset{I}{.Math.}}{\tilde{x}}_i\left({\mathrm{mT}}_r\right){\left(p_k\right)}^i,& \left(47\right)\end{array}$
where {tilde over (x)}.sub.i(mT.sub.r) is the representation of the sum of x(mT.sub.s) and c([(kT.sub.r−mT.sub.s)/T.sub.s,i]). Equation (7) can be implemented by the Farrow structure.

(104) Equation (8) is:

(105) $\begin{array}{cc}\begin{array}{c}y\left(k{T}_r\right)=\underset{m=-\infty }{\overset{\infty }{.Math.}}\underset{i=0}{\overset{I}{.Math.}}x\left(m{T}_s\right)c\left(\left[\frac{k{T}_r-m{T}_s}{T_r}\right],i\right){\left(-\frac{{\mathrm{mT}}_s}{T_r}-\left[-\frac{{\mathrm{mT}}_s}{T_r}\right]-0\mathrm{.5}\right)}^i\\ =\underset{m=-\infty }{\overset{\infty }{.Math.}}\underset{i=0}{\overset{I}{.Math.}}x\left(m{T}_s\right)c\left(\left[\frac{k{T}_r-m{T}_s}{T_r}\right],i\right){\left(p_m\right)}^i.\end{array}& \left(48\right)\end{array}$

(106) As p.sub.m is associated with m, (8) can be simplified to

(107) $\begin{array}{cc}y\left(k{T}_r\right)=\underset{m=-\infty }{\overset{\infty }{.Math.}}\underset{i=0}{\overset{I}{.Math.}}{\hat{x}}_i\left(m{T}_s\right)c\left(\left[\frac{k{T}_r-m{T}_s}{T_r}\right],i\right),& \left(49\right)\end{array}$
where {tilde over (x)}.sub.i(mT.sub.s) is the representation of the sum of x(mT.sub.s) and p.sub.m.sup.i. Equation (9) can be implemented by the Transposed Farrow structure. Equations (6) and (8) can both be denoted as the resampling model based on Farrow structure.

(108) Next, the resampling algorithm based on window function can be used by a broadband frequency measurement device to realize different sampling rate conversions, and provide conditions for the measurement of different frequency components. In addition, the proposed algorithm can flexibly adjust the cutoff frequency of the filter and the window function parameter to prevent interference components in the signals.

IV. Simulation and Application

(109) A. Comparison with Existing Optimal Design Methods without Considering Frequency Aliasing

(110) Because the characteristics of the filter reflect the applicability of resampling, the VFD filter is designed in this example to reflect the design effect of the presented algorithm. Additionally, the presented algorithm is compared with other algorithms under the same conditions. One such algorithm is called WLS and another is called CG. The parameters are set as N=67, I=7, α=0.9π, β=10.5, ω.sub.c=π and δ=0.05. To evaluate the accuracy of the designed VFD filter, some evaluation criteria need to be defined.

(111) (1) The normalized error:

(112) $\begin{array}{cc}{\mathrm{.Math.}}_2={\left(\frac{{\int }_0^{\alpha \pi }{\int }_{-0.5}^{0.5}{.Math.H_p\left(\omega ,p\right)-H_d\left(\omega ,p\right).Math.}^{2}dpd\omega }{{\int }_0^{\alpha \pi }{\int }_{-0.5}^{0.5}{.Math.H_d\left(\omega ,p\right).Math.}^{2}dpd\omega }\right)}^{0.5}\times 100\%,& \left(50\right)\end{array}$
ε.sub.2 reflects the overall magnitude error of the filter passband.

(113) (2) The maximum magnitude error:
ε.sub.∞=20 log.sub.10(max{|H.sub.p(ω,p)−H.sub.d(ω,p)|,ω∈[0,απ],p∈[−0.5,0.5]}),  (51)
ε.sub.∞ reflects the maximum magnitude error of the filter passband.

(114) (3) The maximum group-delay error:
ε.sub.g=max{|τ(ω,p)−p|,ω∈[0,απ],p∈[−0.5,0.5]},  (52)
where p is the group delay of the desired VFD filter and ε.sub.g reflects the maximum group delay error of the filter passband. Simulation tests based on the evaluation criteria are performed, and the design results of three methods are presented in Table II.

(115) TABLE-US-00002 TABLE II COMPARISON BETWEEN ALGORITHMS Method ε.sub.2 (%) ε.sub.∞ (dB) ε.sub.g WLS 2.39 × 10.sup.−4 −99.20 8.90 × 10.sup.−4 CG 2.52 × 10.sup.−4 −101.54 1.36 × 10.sup.−3 Proposed 2.25 × 10.sup.−4 −97.00 2.15 × 10.sup.−3

(116) As can be seen from Table II, compared with other methods, the overall magnitude error ε.sub.2 of the presented method is smaller, but the maximum errors ε.sub.∞ and ε.sub.g are larger. FIG. 4 shows the magnitude error of an exemplary VFD filter. FIG. 5 shows the group delay error of an exemplary VFD filter.

(117) It is evident from FIGS. 4 and 5 that the designed filter characteristics are good in the low and middle-frequency bands, and the maximum magnitude and group delay error both appear at frequency ω=0.9π. Because the maximum magnitude error of the designed filter is −97 dB, the error is close to −100 dB and acceptable. Thus, the method proposed in this paper can be used to design a filter with good performance.

(118) When the designed filter is used in resampling model, it can only be used in the case in which T.sub.s=T.sub.r. To reflect the filter characteristics of the resampling model for other T.sub.r, take T.sub.r=2×T.sub.s and according to (12), the filter order becomes 135. FIG. 6 shows the magnitude responses of the corresponding exemplary VFD filter. In FIG. 6, T.sub.r=2×T.sub.s.

(119) It can be seen from FIG. 6 that the passband range of the filter is halved to [0,0.45π]. However, the magnitude is −6.021 dB when ω=0.5π, and the stopband cutoff frequency ω.sub.s is around 0.55π; therefore, the filter cannot effectively suppress signals with frequency around 0.57. It can be concluded that resampling by the designed filter in this example cannot satisfy the Nyquist theorem, because the cutoff frequency ωc is slightly too large.

(120) B. Comparison with Existing Optimal Design Methods Considering Frequency Aliasing

(121) Because resampling in example A cannot satisfy the sampling theorem due to the high cutoff frequency, the passband and cutoff frequency of the filter are adjusted in this example, and the VFD filter is redesigned by three algorithms.

(122) In this exemplary simulation, the parameters are set as N=67, I=7, α=0.97, β=15, ω.sub.c=0.8571 and δ=0.05. To reflect the filter stopband attenuation effect, a new evaluation criterion is defined:
ε.sub.s=max{20 log.sub.10(|H(π,p)|,p∈[−0.5,0.5])},  (53)
ε.sub.s reflects the suppression effect of the designed filter on signals near ω=π. The comparison results based on the evaluation criteria are presented in Table III.

(123) TABLE-US-00003 TABLE III COMPARISON BETWEEN ALGORITHMS Method ε.sub.2 (%) ε.sub.∞ (dB) ε.sub.g ε.sub.g (dB) WLS —* —* —* —* CG 6.50 × 10−6 −124.00 3.30 × 10−5 −13.31 Proposed 6.14 × 10−6 −124.69 1.78 × 10−6 −153.61  *Fails to converge

(124) The simulation results demonstrate that the method based on WLS cannot design the required filter. The reason for this is that when the passband is narrow, the error between the designed filter and the ideal filter is large, resulting in the failure to converge.

(125) By comparing the four evaluation criteria, it can be found that the characteristics of the VFD filter designed by the proposed method are obviously better than those of the filter designed by CG method. ε.sub.s reveals that CG method cannot effectively suppress signals around ω=π. In comparison, the proposed method can effectively control the stopband of the filter and adjust it according to the requirements. FIG. 7 shows the magnitude error of the designed exemplary VFD filter. FIG. 8 shows the group delay error of the designed exemplary VFD filter.

(126) As can be seen from FIGS. 7 and 8, the VFD filter designed by the presented method has a stable error in the passband, and the group delay error is maintained in the order of 10.sup.−6. The maximum magnitude and group delay error both appear at frequency ω=0.7π, and the maximum magnitude error of the designed filter is −124.69 dB, which demonstrates that the magnitude error in passband is very small.

(127) FIG. 9 shows the magnitude responses of the exemplary VFD filter when T.sub.r=2×T.sub.s. Like in example A, this reflects the filter characteristics of the resampling model for other T.sub.r.

(128) It can be seen from FIG. 9 that the passband range of the filter is halved to [0,0.35π] and the stopband cutoff frequency is around ω=0.57c where the magnitude is about −159.6 dB. Thus, in this example, the designed filter can filter out signals above 0.5π, and resampling can satisfy the Nyquist theorem when T.sub.r=2×T.sub.s.

(129) Therefore, it can be concluded that the proposed method is more adaptable than the traditional optimal method, and the method can flexibly adjust Kaiser window parameter p and cutoff frequency ω.sub.c to meet the needs of practical applications.

(130) C. Signal Processing with Resampling Model

(131) To verify the effect of the resampling algorithm when applied to the broadband frequency measurement device, this example applies the designed subfilter coefficients in example B to resampling model (8), uses (8) to resample the signal at different sampling rates, and then analyzes the resampling results.

(132) The original signal with harmonic and inter-harmonic components is considered, let
x(t)=√{square root over (2)}[57.73 cos(100πt)+5.773[cos(200πt)+cos(6190πt)]].  (54)

(133) For the original signal x(t), sampling is performed at the sampling rate of f.sub.s=50 kHz. The sampling rate is then converted to f.sub.r, the resampling values are output by the resampling model, and the amplitudes of different frequency signals are finally analyzed by Fast Fourier Transform. The results are exhibited in Table IV.

(134) TABLE-US-00004 TABLE IV RESULTS OF THE RESAMPLING ALGORITHM AT DIFFERENT SAMPLING RATES f.sub.r Method 50 Hz (V) 100 Hz (V) 2905 Hz (V) 3095 Hz (V) 50 kHz CG 57.73 5.773 0 5.773 Proposed 57.73 5.773 0 5.773 10 kHz CG 57.73 5.773 0 5.773 Proposed 57.73 5.773 0 5.773  6 kHz CG 57.73 5.773 0.059 / Proposed 57.73 5.773 1.6 × 10.sup.−10 /

(135) As presented in Table IV, when f.sub.r=6 kHz, signals of −2905 Hz and 2905 Hz exist simultaneously in the design results by CG method, resulting in some interference signals after resampling. Compared with CG method, the presented method can effectively suppress the signals around 3000 Hz and prevent the appearance of interference signals. Additionally, the resampling model does not affect the low-frequency signals far away from f.sub.r/2, and the filter in the model suppresses the signals near f.sub.r/2.

(136) It is also necessary to verify that the resampling has the function of eliminating the deviation of the non-synchronous sampling to achieve synchronous sampling. This example uses the resampling model (8) to sample signals with different delays, and the results are then analyzed.

(137) The original signal with harmonics and delay is considered, let
x(t)=√{square root over (2)}[57.73 cos(100π(t−(33.5+p)/f.sub.s))+5.773[cos(500π(t−(33.5+p)/f.sub.s)),+cos(5000)(t−(33.5+p)/f.sub.s))]]  (55)
where t-33.5-p is considered the group delay of the fractional delay filter. For the original signal x(t), sampling is performed at the sampling rate of f.sub.s=f.sub.r=50 kHz, and the initial phase angles of different frequency signals in the resampling results are analyzed by Fast Fourier Transform. In addition, the indicator p.sub.d is used to reflect the phase angle deviation caused by the lack of compensation for p, and the indicators p.sub.c and p.sub.w respectively represent the phase angle error compensated by the resampling models of CG algorithm and the proposed algorithm. The results are presented in Table V.

(138) TABLE-US-00005 TABLE V SAMPLING RESULTS OF THE RESAMPLING ALGORITHM AT DIFFERENT DELAYS Phase angle p error 50 Hz (°) 250 Hz (°) 2500 Hz (°) −0.5 p.sub.d 0.18 0.9 9 p.sub.c 3.4 × 10.sup.−7 1.4 × 10.sup.−6 −1.4 × 10.sup.−6   p.sub.w −1.1 × 10.sup.−8   −4.8 × 10.sup.−8   4.5 × 10.sup.−8 −0.3 p.sub.d 0.108 0.54 5.4 p.sub.c −2.1 × 10.sup.−8   −8.5 × 10.sup.−8   6.0 × 10.sup.−8 p.sub.w −5.7 × 10.sup.−8   −2.4 × 10.sup.−7   2.7 × 10.sup.−7 0 p.sub.d 0 0 0 p.sub.c .sup. −6.1 × 10.sup.−10   −6.8 × 10.sup.−9   .sup. −6.5 × 10.sup.−10   p.sub.w .sup. −6.1 × 10.sup.−10   .sup. −6.8 × 10.sup.−10   .sup. −6.5 × 10.sup.−10   p.sub.d −0.072 −0.36 −3.6 0.2 p.sub.c 1.7 × 10.sup.−8 6.1 × 10.sup.−8 −5.6 × 10.sup.−8   p.sub.w 4.9 × 10.sup.−8 2.0 × 10.sup.−7 −2.3 × 10.sup.−7   p.sub.d −0.18 −0.9 −9 0.5 p.sub.c −3.4 × 10.sup.−7   −1.4 × 10.sup.−6   1.4 × 10.sup.−6 p.sub.w −9.5 × 10.sup.−9   3.5 × 10.sup.−8 −4.6 × 10.sup.−8  

(139) The results in Table V demonstrate that the resampling model can compensate for different delays p. In contrast, the compensation effects of the two methods are similar.

(140) Therefore, the method proposed in this disclosure has sound results in signal processing, and its function meets the conditions for application in broadband frequency measurement device.

V. Conclusion

(141) This disclosure describes a resampling algorithm based on window function, that can provide conditions for different sampling rates required by existing measurement equipment to measure different component signals. The frequency responses of the filter in resampling model are revealed, and the design conditions of the filter are specified. Then a fractional delay filter design model based on window function is established. According to the model, fractional delay filter coefficients with different group delay responses are designed, and a filter coefficients matrix is established. The polynomial form of the filter coefficients matrix is established, and subfilter coefficients are solved by least square method. As verified by simulation tests, the proposed method can flexibly adjust the cutoff frequency of the filter and the window function parameter to prevent frequency aliasing. When the filter passband requirement is small, the presented method has a better effect than the existing methods. In addition, the presented algorithm can meet the requirements for sampling rate conversion and delay compensation, and can be applied to broadband frequency measurement devices.

(142) Technical Implementation of the Server

(143) FIG. 10 illustrates exemplary hardware components of a server. A computer system 1000, or other computer systems similarly configured, may include and execute one or more subsystem components to perform functions described herein, including the steps of various flow processes described above. Likewise, a mobile device, a cell phone, a smartphone, a laptop, a desktop, a notebook, a tablet, a wearable device, a server, etc., which includes some of the same components of the computer system 1000, may run an application (or software) and perform the steps and functionalities described above. Computer system 1000 may connect to a network 1014, e.g., Internet, or other network, to receive inquiries, obtain data, and transmit information and incentives as described above.

(144) The computer system 1000 typically includes a memory 1002, a secondary storage device 1004, and a processor 1006. The computer system 1000 may also include a plurality of processors 1006 and be configured as a plurality of, e.g., bladed servers, or other known server configurations. The computer system 1000 may also include a network connection device 1008, a display device 1010, and an input device 1012.

(145) The memory 1002 may include RAM or similar types of memory, and it may store one or more applications for execution by processor 1006. Secondary storage device 1004 may include a hard disk drive, floppy disk drive, CD-ROM drive, or other types of non-volatile data storage. Processor 1006 executes the application(s), such as those described herein, which are stored in memory 1002 or secondary storage 1004, or received from the Internet or other network 1014. The processing by processor 1006 may be implemented in software, such as software modules, for execution by computers or other machines. These applications preferably include instructions executable to perform the system and subsystem component functions and methods described above and illustrated in the FIGS. herein. The applications preferably provide graphical user interfaces (GUIs) through which users may view and interact with subsystem components.

(146) The computer system 1000 may store one or more database structures in the secondary storage 1004, for example, for storing and maintaining the information necessary to perform the above-described functions. Alternatively, such information may be in storage devices separate from these components.

(147) Also, as noted, processor 1006 may execute one or more software applications to provide the functions described in this specification, specifically to execute and perform the steps and functions in the process flows described above. Such processes may be implemented in software, such as software modules, for execution by computers or other machines. The GUIs may be formatted, for example, as web pages in HyperText Markup Language (HTML), Extensible Markup Language (XML) or in any other suitable form for presentation on a display device depending upon applications used by users to interact with the computer system 1000.

(148) The input device 1012 may include any device for entering information into the computer system 1000, such as a touch-screen, keyboard, mouse, cursor-control device, microphone, digital camera, video recorder or camcorder. The input and output device 1012 may be used to enter information into GUIs during performance of the methods described above. The display device 1010 may include any type of device for presenting visual information such as, for example, a computer monitor or flat-screen display (or mobile device screen). The display device 1010 may display the GUIs and/or output from sub-system components (or software).

(149) Examples of the computer system 1000 include dedicated server computers, such as bladed servers, personal computers, laptop computers, notebook computers, palm top computers, network computers, mobile devices, or any processor-controlled device capable of executing a web browser or other type of application for interacting with the system.

(150) Although only one computer system 1000 is shown in detail, system 1000 may use multiple computer systems or servers as necessary or desired to support the users and may also use back-up or redundant servers to prevent network downtime in the event of a failure of a particular server. In addition, although computer system 1000 is depicted with various components, one skilled in the art will appreciate that the system can contain additional or different components. In addition, although aspects of an implementation consistent with the above are described as being stored in a memory, one skilled in the art will appreciate that these aspects can also be stored on or read from other types of computer program products or computer-readable media, such as secondary storage devices, including hard disks, floppy disks, or CD-ROM; or other forms of RAM or ROM. The computer-readable media may include instructions for controlling the computer system 1000, to perform a particular method, such as methods described above.

(151) The present disclosure is not to be limited in terms of the particular embodiments described in this application, which are intended as illustrations of various aspects. Many modifications and variations can be made without departing from its spirit and scope, as may be apparent. Functionally equivalent methods and apparatuses within the scope of the disclosure, in addition to those enumerated herein, may be apparent from the foregoing representative descriptions. Such modifications and variations are intended to fall within the scope of the appended representative claims. The present disclosure is to be limited only by the terms of the appended representative claims, along with the full scope of equivalents to which such representative claims are entitled. It is also to be understood that the terminology used herein is for the purpose of describing particular embodiments only, and is not intended to be limiting.