I am trying to implement the following 4-Tap Polyphase window shown in figure below and used in : polyphase filters- radio astronomy: lecture #8 by Prof. Dale E. Gary on Matlab. I have tried splitting the window in time domain after multiplication with Sinc function but has not achieved the sidelobe level and mainlobe width shown in green graph- figure on right. I am not sure what I might be missing here.
I assume that the Polyphase filter shown in the paper is used as a decimator such that the output rate is an integer fraction of the input rate (and this would make sense in a correlation operation since the output frequency required is less, due to the proper filtering of the correlator.) Making a polyphase filter implementation is quite easy; given the desired coefficients for a simple FIR filter, you distribute those same coefficients in "row to column" format into the separate polyphase FIR components as explained in the following example:
Assume a FIR filter with 8 taps for simplicity of explanation, with coefficients as follows:
h[0], h[1], h[2], h[3], h[4], h[5], h[6], h[7], ...
and we use that to create a 4 element polyphase each with 2 taps:
Poly FIR1: h[0], h[4]
Poly FIR2: h[1], h[5]
Poly FIR3: h[2], h[6]
Poly FIR4: h[3], h[7]
From the above it is hopefully clearer in what I meant by "row to column" format.
The Polyphase Decimator would be structured as in the figure below:
A Polyphase Interpolator works the same way.The graphics below explain clearly why we do the row to column mapping. This shows an input signal of successive pulses at a higher rate going through an interpolation operation. Normally, without the use of polyphase implementations, we can interpolate a signal by simply inserting zeros, and then following that with a low pass filter to get rid of the higher frequency images that occur due to the zero insert. The best filter design would pass our original spectrum with no distortion and completely eliminate the higher frequency components. It is the coefficients of this FIR filter that are used in the polyphase implementation since it mathematically performs the same result (that said, I suspect the filter chosen in the paper which used the windowing method of filter design is not the optimum choice compared to a least squares design solution).
The technique in the paper may be misnamed (or does not fit the normal use of polyphase filtering for resampling).
Normally when a window is made shorter than the FFTs length (by zero-padding, etc.), the main lobe of the frequency response of the windowing artifact gets wider.
The paper's FFT filter seems to be using the technique of making the window on the data be longer than the FFTs length. How do you fit a longer window into a shorter FFT? By (additively) wrapping the windowed data circularly around the FFT vector. That narrows the main lobe of the frequency response, assuming stationarity across a width longer than the FFT, and thus at the cost of time locality related to the FFT's width. Also note that some frequency sinusoids can cancel themselves out using this method, as the windowed FFT bin width can become narrower than the FFT bin spacing.
The technique in the paper seems to be wrapping the pre-windowed data around the FFT input vector 4 times, additively.
You can use this code to perform tests (for Matlab or Octave). This basically writes down two sinusoids, and analyzes them following the very procedure as in https://casper.berkeley.edu/wiki/The_Polyphase_Filter_Bank_Technique
Keep in mind that the base unit for the frequency is "df" (in the code below); a frequency in the signal, that is an integer multiple of "df", will have no leakage (and the so-called polyphase technique will be of no use). In the code you chose f = (integer)*df for no leakage, or (integer-decimal_point-some_number)*df to have leakage (a real situation). You should find that this technique leaves the actual spectral resolution unchanged, but strongly reduces the leakage; so you will be able to spot a weak f-peak close to a strong f-peak. With frequencies f1 = 17.3*df and f2 = 21.3*df, and the f1/f2 sinusoids amplitudes 3.4/0.8 (just as an example), you will see the effect. In practice, you can hardly tell different peaks from each other, in either way, if they are spaced less than (say) 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)');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
