我假设论文中显示的多相滤波器用作抽取器，使得输出速率是输入速率的整数分数（这在相关运算中是有意义的，因为由于适当的过滤，所需的输出频率较小相关器。实现多相滤波器非常容易；给定简单 FIR 滤波器所需的系数，您可以将这些相同的系数以“行到列”格式分配到单独的多相 FIR 组件中，如下例中所述：

为简化说明，假设 FIR 滤波器有 8 个抽头，其系数如下：

h[0], h[1], h[2], h[3], h[4], h[5], h[6], h[7], ...

我们用它来创建一个 4 元素多相，每个多相有 2 个抽头：

聚 FIR1：h[0]、h[4]

聚 FIR2：h[1]、h[5]

聚 FIR3：h[2]、h[6]

聚 FIR4：h[3]、h[7]

从上面可以看出我所说的“行到列”格式的含义会更清楚。

多相抽取器的结构如下图所示：

多相插值器的工作方式相同，因为我已经有图表清楚地解释了为什么我们进行行到列的映射，我将在下面包含它（以及直观地解释多相实现本身的美！）。这显示了以较高速率通过插值操作的连续脉冲的输入信号。通常，在不使用多相实现的情况下，我们可以通过简单地插入零来对信号进行插值，然后使用低通滤波器来消除由于插入零而出现的高频混叠。最好的滤波器设计将通过我们的原始频谱而不会失真，并完全消除高频分量。

论文中的技术可能命名错误（或不适合正常使用多相滤波进行重采样）。

通常，当窗口短于 FFT 长度（通过零填充等）时，窗口伪影的频率响应的主瓣变得更宽。

该论文的 FFT 滤波器似乎使用了使数据窗口长于 FFT 长度的技术。如何将较长的窗口放入较短的 FFT 中？通过（附加地）将窗口数据循环包裹在 FFT 向量周围。这会缩小频率响应的主瓣，假设在比 FFT 更长的宽度上具有平稳性，因此以与 FFT 宽度相关的时间局部性为代价。另请注意，某些频率正弦曲线可以使用此方法自行抵消，因为窗口化 FFT 箱宽度可能变得比 FFT 箱间距更窄。

论文中的技术似乎是将预加窗数据围绕 FFT 输入向量附加 4 次。

您可以使用此代码执行测试（对于 Matlab 或 Octave）。这基本上写下两个正弦曲线，并按照https://casper.berkeley.edu/wiki/The_Polyphase_Filter_Bank_Technique中的程序分析它们

请记住，频率的基本单位是“df”（在下面的代码中）；信号中的频率，即“df”的整数倍，将没有泄漏（所谓的多相技术将没有用处）。在代码中，您选择 f = (integer)*df 表示没有泄漏，或 (integer-decimal_point-some_number)*df 表示有泄漏（真实情况）。您应该会发现，这种技术使实际光谱分辨率保持不变，但大大减少了泄漏；因此您将能够发现接近强 f 峰值的弱 f 峰值。频率 f1 = 17.3*df 和 f2 = 21.3*df，f1/f2 正弦波幅度为 3.4/0.8（仅作为示例），您将看到效果。在实践中，如果它们的间距小于（例如）1.5*df，无论哪种方式，您都几乎无法区分不同的峰值。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all
close all

N = 300;%number of points for the signal

T = 1;%1.9531e-3;%sampling time

t = 0:T:(N-1)*T;%time vector

df = 1/(N*T);%base frequency

f1 = 17.3*df;%example freq 1

f2 = 21.3*df;%example freq 2


x = 3.4*sin(2*pi*f1*t) + 0.8*sin(2*pi*f2*t);%the signal (two close freq)

for i=1:N,
  x(i) = x(i) +2*(rand-0.5);%noise added

endfor;

%plot(t,x);

%here the same signal, but 4 times long as much, in time

t_e = 0:4*N-1;
x_e = 3.4*sin(2*pi*f1*t_e) + 0.8*sin(2*pi*f2*t_e);
for i=1:4*N,
  x_e(i) = x_e(i) +2*(rand-0.5);%noise added

endfor;

w = blackman(4*N);
for i=1:4*N,
  x_e(i) = x_e(i)*w(i);%windowing

endfor;

%plot(t_e,x_e);

%here the 4-fold wrapping:
for i=1:(4*N)/4,
  x1(i) = x_e(i);
endfor;
for i=1:(4*N)/4,
  x2(i) = x_e(i +(4*N)/4);
endfor;
for i=1:(4*N)/4,
  x3(i) = x_e(i +2*(4*N)/4);
endfor;
for i=1:(4*N)/4,
  x4(i) = x_e(i +3*(4*N)/4);
endfor;

xp = x1 +x2 +x3 +x4;%time wrapping


%plot(t,x4);

%plot(t, xp);

fstep = 1/(N);    %FFT frequency step
f = zeros(1,round(N/2));
for i = 1:round(N/2),%frequency vector
  f(i) = i*fstep-fstep;
end


Y = 2*fft(x)/N;% FFT transform of the short signal
Y = Y(1:round(N/2));%first half of the spectrum selected (the only useful, as the signal is real)
Y2 = abs(Y).^2;%squared amplitude (power)
%plot(f,real(Y),'g'); hold on;%parte reale della FFT
%plot(f,imag(Y),'r'); hold on;%parte immaginaria della FFT
%plot(f,Y2); hold on; xlabel('Frequenza (Hz)');%valore assoluto (alquadrato) della FFT
plot(f,Y2,'r-o'); hold on; xlabel('Frequenza (Hz)');


Yp = 2*fft(xp)/(N);%FFT transform of the long, wrapped, signal
Yp = Yp(1:round(N/2));%selezione della prima metà  dello spettro (unicoutile, se P è reale)
Y2p = abs(Yp).^2;%ampiezza reale quadratica dello spettro (potenza)
%plot(f,real(Y),'g'); hold on;%parte reale della FFT
%plot(f,imag(Y),'r'); hold on;%parte immaginaria della FFT
%plot(f,Y2); hold on; xlabel('Frequenza (Hz)');%valore assoluto (al quadrato) della FFT
plot(f,Y2p,'g-o'); hold on; xlabel('Frequenza (Hz)');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%