6 tap 32 phase polyphase interpolation

[trevlac]

I finally figured out what this means. The BT878 has a FIR filter with 192 coefficients. That's a bunch. That is also why it does a mean resize. I think it has a cutoff of about 5MHz, but that's a guess. It is definately one of those 'sin c' graphs.

6 * 32 = 192

When you do polyphase you divide the taps by the phases so you get 32 phases with 6 taps each. Then you can decimate 1st and apply the filter to every x pixels, up to 32.

IE: Resize by 2/3 from 720 to 480.


Upsize by 2 putting a zero pixel between each real pixel = 1440 pixels
Seperate the 1440 into 1440/3 lines of 480 each. Throw out the last line
Seperate your 192 tap filter into 3 64 tap filters, Throw out the last
Convolute the 2 filters with the 2 lines and sum the results


You now have decimated by 2/3 by using a 192 tap FIR filter.

I'd say this would be quite slow, but an interesting experiment. I'm not sure how to get the coefficients for such a filter. I've found tools that do up to 80.


PS: Divide/seperate out means if you have pixles 1,2,3,4,5,6,7,8,9 seperate/divide them into groups like this: 1,4,7 : 2,5,8 : 3,6,9 to 'divide' by 3.

[Arachnotron]

Is it not by simply calculating the value of sinc(x) for all 6x32 points? (In fact, since sinc is symmetrical, you should be able to get by using only 96 values)

I have been messing around a bit too. Trying to write a dummy filter working only on a raw 8 bit grayscale file using a Kaiser Window, which is one of the solutions for the aliassing problems you get using a truncated sinc function. (sinc runs to + and - infinite, and you have to stop somewhere. This causes aliassing effects. Lancoz and Kaiser are adjustments to a pure sinc function to minimize this effect)

[trevlac]

Originally posted by Arachnotron
Is it not by simply calculating the value of sinc(x) for all 6x32 points? (In fact, since sinc is symmetrical, you should be able to get by using only 96 values)


That's easy for you to say. :D

If you use different windows you get very different results close to the cutoff point.

See this calculator:

http://www.ittc.ku.edu/~rvc/projects/firfilter/windowed/

And you are correct on the 96 value thing. The above calculator only does 40 values. You then double it. To get the 80 they talk about.

I am going to try to write an avisynth resize filter to use maybe 64 taps total (32 values). Any suggestions for a cutoff near 50% of the sampling frequency? I can't say I really understand this.

Trev

[Arachnotron]

That's easy for you to say

My understanding:

As I understand it, you have a window stretching from -z to + Z on either side of the pixel you are looking at. You choose the function you want to use - lancoz, Kaiser or simply a truncated sinc - and calculate the value of this function for the 6 points between the samples -Z and -Z+1, -Z+1 and -Z+2 .... Z-1 and Z. This gives you an 6x32=192 wide array of function values. (or actually you use only half because of the symetric function) Then you move this array down the line, and summarize the contribution made by each old sample between -Z and +Z for each new sample between the old samples.

At this point you probably need to do some normalization.

At the edges you have a problem, since there are no values on the left or the right. So you can either take Z fake 50% grey pixels on the left / right OR put in the average luma value of the last Z pixels you have. There are probably other methods for that too, but I haven't fount them yet :)

All this Kernel stuff is just a re-write of this principle for speed optimization.

From the litt, I get the impression Kaiser is the best, but not so fast because a good implementation uses floats.

A kaiser function has two coefficients, Tau and Alpha. Tau is actually just Z, the width of your window. In combination with Alpha this will define the width of your stopband. What the exact relation is between alpha and stopband width I still don't understand. I have a link somewhere of an article that discusses the different reconstruction filters. Will post it tomorrow.

[trevlac]

Originally posted by Arachnotron
As I understand it, you have a window stretching from -z to + Z on either side of the pixel you are looking at. You choose the function you want to use - lancoz, Kaiser or simply a truncated sinc - and calculate the value of this function for the 6 points between the samples -Z and -Z+1, -Z+1 and -Z+2 .... Z-1 and Z. This gives you an 6x32=192 wide array of function values. (or actually you use only half because of the symetric function) Then you move this array down the line, and summarize the contribution made by each old sample between -Z and +Z for each new sample between the old samples.


I don't completely follow. Should it not be 32 points between the samples? And are you saying each point gets it's own sinc filter?

For example, for a 5 tap 2 phase:

2 sinc filters
0.5,-0.5,1,-0.5,0.5 and 0.5,-0.5,1,-0.5,0.5

or

1 sinc filter
-0.08, -0.11, 0.0, 0.21, 0.41, 0.41 , 0.21, 0.0, -0.11, -0.08

Which get changed to :
-0.08,0.0, 0.41, 0.21, -0.11
-0.11,0.21,0.41, 0.0,-0.08

Maybe they are the same, but I have a bad example.

Using more taps does seem to allow a sharper cutoff. The edges are indead a problem. Possibly up sizing a lot, to reduce the effect of the grey on the edges is what they do. For 1/2 go up by 16. That way you only ever have 12 taps to go over the edge. This would effectively use only 2*12 of the 192 coefficients because the others would be on zero lines. I guess. :)

I really have to try this to get it down.

[Arachnotron]

I still get confused about the term tap and phase. Perhaps I should just drop it :)

I am simply trying to get the algorithm working, completely inefficient, and worry about optimization later.

suppose I want to interpolate a new point between two samples. I select the number of original samples on either side I want to contribute to this new point. I thought these were my phases, but apparently not? Suppose Z old samplepoints on either side, so 2*Z total.

Anyhow, I select my function - Kaiser in this case-
and take pixel -Z on the left of the new sample. I calculate the value of the Kaiser function for sample -Z at the new point, multiply it with the luma value of point -Z and drop it in the new sample.
Next I move to point -Z+1. Now calculate the Kaiser function of point -Z+1 at the same new samplepoint, multiply with luma of sample -Z+1 and add this to the luma value of the new point. Repeat until point +Z.
The final luma value of the new point is now made up of contributions of all -Z to Z samples.

Now move over to the next new sample point you want to know about and repeat the procedure. Since the new points are equidistant, you will only need to calculate the Kaiser function for certain points between the old samples and dump them in a lookup table. I thought these new points were the taps?

So in fact you only need to calculate the Kaiser function for Z*(number of new points between samples). And normalize them, so the sum of all the values in this window is 1.

I thought the values in this lookup table are the coefficients you use?

I try to implement the Kaiser function from this (http://citeseer.nj.nec.com/301662.html) page, see pdf link at top right corner, "Mastering Windows: Improving Reconstruction". Will take me a while; perhaps this weekend.

[edit] this is a load if bullsh***t. :( For info on filters, see
this (http://virtualdub.everwicked.com/index.php?act=ST&f=2&t=5729&) thread over at everwicked

[trevlac]

With phaeron's explaination and some info Arachnotron provided, I'd like to give this 1 more try. :D In the effort to learn just a bit more.

General process then an example:

1) Determine the decimation factor as 1/scale ie: 2/3 scale = 3/2.

2) Plug (phases - 1) 'zero' pixels between your source pixels

3) Keep the 1st source pixel as your 1st dest pixel. Then count across this line and find the next (phase * decimation factor) pixel. ie: 8 phase * 3/2 = 12th pixel.

4) The 12th pixel should fall as one of 8 between 2 real pixels. Find the filter you use for this phase. You should have 8 filters to choose from for an 8 phase process.

5) Filter using the surounding Xtap real pixels. ie: if you have a 2 tap filter, use the surounding 2 real pixels. The 1st source pixel should also be filtered the same way.


Example of 8 phase 2 tap decimation by 2/3
Go from 6 pixels to 4.

1) Decimation factor is 3/2.
source is ABCDEF, destination is abcd

2) A2345678B2345678C2345678D2345678E2345678F2345678

3/4) A1 B5 D1 E5
assume filters are
.27,.02
.56, .06
.36, 0
.12,.04
.16,.01
.09,.03
.08,.01
.08,.02

So we want the use .27,.02 on A1 and D1 and .16,.01 on B5 and E5.

5)
A1 = convolution of filter(.27,.02) and A + (the pixel before A)
B5 = convolution of filter(.16,.01) and A + B
D1 = convolution of filter(.27,.02) and C + D
E5 = convolution of filter(.16,.01) and D + E


----------
Some things to question:
1) Are the filters their own, or are they a larger filter broken into 8 parts?


I may be babbeling, but I'm really trying to figure this out.
:confused:

[Wilbert]

A1 B5 D1 E5
I don't understand these numbers :)

of 8 phase 2 tap decimation by 2/3
As I understand it, 8 phase filter indicates the maximal number of phases or sub-filters. Thus, interpolating by a factor of 8 is the maximal.

I'm not really sure whether the above is true, will ask that in the thread in the vdub forum.

In this case you start with interpolating by a factor of 2, which means that you only use 2 phases or sub-filters:

A B C D E F => A0B0C0D0E0F

two sub-filters: (h0,h2), (h1,h3)

(thus the corresponding 'total' FIR filter is 4 taps with this interpolation)

result:
h0*0 + h2*A = x0
h1*0 + h3*A = x1

h0*A + h2*B = x2
h1*A + h3*B = x3

h0*B + h2*C = x4
h1*B + h3*C = x5

h0*C + h2*D = x6
h1*C + h3*D = x7

h0*D + h2*E = x8
h1*D + h3*E = x9

h0*E + h2*F = x10
h1*E + h3*F = x11

That is 12 new pixels. Now you decimate by a factor of 3, keeping each third sample:

h0*A + h2*B = a
h1*B + h3*C = b
h0*D + h2*E = c
h1*E + h3*F = d

Looks a bit similar as what you have, but I don't understand the numbers in filter(.27,.02) ... Or are they computed form the sinc function?

Btw, Avery said this in a thread (bt8x8 chip):
"the decimation filters aren't specified but my guess is that they're simply a set of Gaussian filters: [1], [1 1]/2, [1 2 1]/4, [1 3 3 1]/8. This would be consistent with the vertical decimation filters."

I guess, this is not equivalent with keeping each third sample. If you decimate by a factor of 3, you use this filter [1 2 1]/4. Result is:

(x0+2*x1+x2)/4
(x3+2*x4+x5)/4
(x6+2*x7+x8)/4
(x9+2*x10+x11)/4

??? Could be totally wrong though?

[trevlac]

Read this: http://www.intersil.com/data/AN/an9661.pdf


Polyphase Filter Design

The coefficients for the interpolation phases are generated
by designing the filter at the interpolated rate (32x the input
sample rate ..... the prototype filter would have (32 phases x 8 taps)
256 taps. These taps are divided into the 32 filter phases by
taking every 32nd sample and storing the result as
coefficient sets 0 through 31. The first coefficient set would
be C0, C32, C64, C96, C128, C160, C193, and C224. The
last coefficient set would be C31, C63, C95, C127, C159,
C192, C223, and C255.

That sounds like they design a big 256 coefficient sinc filter and break it up.

Also, in secion "3.4.2 Can you give me an example of a FIR interpolator?" : here http://www.dspguru.com/info/faqs/multrate/interp.htm

They clearly discuss a 'longer' fir that is broken up as I've described.


Finally, below is a book on the topic. I'd buy it if I wasn't affaid that I would not understand a word. :)

http://www.amazon.com/exec/obidos/tg/detail/-/0136057187/ref=pd_bxgy_text_1/103-5930648-4836658?v=glance&s=books&st=*

[Wilbert]

I have seen your first link, and it doesn't talk about sinc interpolation. I will look for the book in our library.

In the mean time, I still hope for an explanation how you calculated the numbers

3/4) A1 B5 D1 E5
assume filters are
.27,.02
.56, .06
.36, 0
.12,.04
.16,.01
.09,.03
.08,.01
.08,.02

If I use equation 11 from this page
http://ccrma-www.stanford.edu/~jos/resample/Appendix_Periodic_Sinc_Interpolation.html

I get x0, 0.97*x0, 0.139*x1, 0.065*x2, 0.042*x3, 0.031*x4, ... for t=0/8, 1/8, 2/8, etc. It doesn't look anywhere near your coefficients.

@Arachnotron,
[edit] this is a load if bullsh***t. For info on filters, see
this thread over at everwicked
Thanks for that thread. But note that Vdub doesn't use sinc interpolation.

Do they mention sinc interpolation in the specs of bt8x8?

[Arachnotron]

@Wilbert
Thanks for that thread. But note that Vdub doesn't use sinc interpolation. Do they mention sinc interpolation in the specs of bt8x8?

I was just referring to my attempt to explain an interpolation filters coefficients, not to whether or not the BT878 uses a truncated sinc function.

The only description I could find about the BT878 (not from the datasheet by the way, but an application note) mentiones the following:

The algorithm uses polyphase interpolation to generate the scaled pixels from the 4xFsc sampled pixels. The desired sampled pixel location is determined to within 1/32 of a 4xFsc-sample period (the Bt835 has 1/64th resolution for luma). These values, between 1 and 32, are used to determine the coefficients used by the 6-tap weighted sum interpolator to calculate the output pixel luminance value. For chroma data, the same 1 of 32 samples are used, but fewer coefficients are used by a 2-tap interpolator.

More than this I could not find. (apart from the formulas in fig. 2.4 of the fusion datasheet, but those don't tell much either)

[trevlac]

Originally posted by Wilbert

In the mean time, I still hope for an explanation how you calculated the numbers

3/4) A1 B5 D1 E5
assume filters are
.27,.02
.56, .06
.36, 0
.12,.04
.16,.01
.09,.03
.08,.01
.08,.02



I'm terribly sorry. The coefficients in my example were given in a rush, and have some problems. They were not the point of my post, but just an example. I think i used this (http://www.ittc.ku.edu/~rvc/projects/firfilter/windowed/) to generate them. Here (http://www.poynton.com/Poynton-dsp.html) is a collection of links to programs that generate FIR filters.

I don't think the BT uses a sinc filter, but an approximation. I think the BT uses 192 coefficients. I see it as marginally better than VDub. This may be due to the choice of coefficients. I can not follow / match up what Avery said with how it appears that polyphase interpolation works.

The point of my post was (for 32 phase 6 tap interpolation); you take a big 192 coefficient filter, break it up into 32 groups with 6 coeffieicents each. Then when you resize, you zero stuff 31 points between each of your original pixels. One originial pixel + these 31 points = the phase/groups of your filter. To decimate back down you count across your new line finding every 1/scale rate * 32 point. For example if you scale by 2/3 (720->480) you find every 3/2*32=48th point. This point tells you which of the 32 6 tap groups to use to calculate the new pixel. The new pixel is convoluted from the original proceding pixels and the given 6 tap group (just like normal convolution).

I will try to do another example. I have change my notation a little to make it more clear. I will also use a 'normal looking' set of coefficients, but I am just making them up. No calculation.

----------------------

8 phase 2 tap convolution of 9 pixels to 6 pixels

Starting Pixels:
A B C D E F G H I J

Please note, we will start with B because a 2 tap filter requires pixels before the pixel we are working on. A is before B, so B works. I don't know what is before A :)


Our big 8*2 coefficient filter
0 ,0.0001,-0.0009,0.0048 ,-0.0187,0.0576 ,-0.16 ,0.617,
0.617,-0.16 ,0.0576 ,-0.0187,0.0048 ,-0.0009,0.0001,0

I found I could not come up with one that made sense, so I used this (http://maxflat.ec.ss.titech.ac.jp/PlotMFnew.htm). N(order)=15, K(?)=1, d(?)=0/1. I then rounded the results to 4 decimals.

Step #1 - Break up our big filter into 8 groups
G1 = 0,0.617
G2 = .0001,-0.16
G3 = -.0009,.0576
G4 = .0048,-0187
G5 = -.0187,.0048
G6 = .0576,-.0009
G7 = -.16,.0001
G8 = .617,0

Step #2 - Zero stuff the pixel line. 8 phase means 7 zeros between each
Here our original line will be denoted as B1, C1, D1 etc.
The new 'zero' pixels are b2,b3,b4,b5,b6,b7,b8,c2,c3, etc.

B1 b2 b3 b4 b5 b6 b7 b8 C1 c2 c3 c4 c5 c6 c7 c8 D1 d2 d3 ..... J1

Step #3 - Find the points we are interested in.
9/6 = 3/2. 3/2*8 phase = Every 12th point.

B1, c5, E1, f5, H1, i5 are the output points, so we use filters G1 and G5 to get them.

Step #4 - Convolution (*)

Our New Output Pixels (U V W X Y Z)

U (B1) = G1 (*) {B,A}
V (c5) = G5 (*) {D,C}
W (E1) = G1 (*) {E,D}
X (f5) = G5 (*) {G,F}
Y (H1) = G1 (*) {H,G}
Z (i5) = G5 (*) {J,I}


Note Next to the output pixels, I have identified the 'point' they came from. I used this to determine the filter used, and also the input pixels used. Honestly, I'm not sure on the input pixels to the convolution. For convolution I think you flip the pixel order (or the filter). So hence my notation of G1*{B,A} for the 1st output pixel. However, I'm guessing that you center the filter on the new pixel and take input pixels from both sides.

PS
Since you are both still interested in this, I will try to write a program to simulate and test different sets of input coefficients.