DCT-IV, DST-IV and DCT-I, DST-I orthonormalization support in scipy.fftpack

[zywina]

Implementation of norm='ortho' option for type 1 dct and dst.

I would appreciate if someone with access to Matlab could produce an equivalent to gendata.m and test.mpz in the test directory to improve the unit tests.

edit: added DCT-IV and DST-IV, description below

[ilayn]

@zywina Thanks for the PR. Is there any example that the tests are written in such a way that the result is known beforehand? Relying on matlab is ok of course but still makes it very difficult to debug if something else breaks this function and we would need matlab once again to track where it goes wrong.

Also I am curious if contributing to FFTPACK directly would be a better way to add this feature since in case we update our FFTPACK version this has to be done again.

[zywina]

I've checked the results against the textbook definitions. Below is some sample code. All of the Fourier related transforms can be reduced to a matrix multiplication. All that needs to be done is to construct the appropriate matrix and verify that the transform produces the same outputs. This sort of approach may be a better long term option for unit tests as it could be applied to other transforms and use randomized inputs.

I don't believe contributing to FFTPACK is viable as the project is inactive. The last release was in 2011 and the version used in scipy is much older than that, perhaps from the 90s.

from scipy import *
from scipy.fftpack import *

def dct1_matrix(N):
    M = zeros((N,N))
    for i in range(N):
        for j in range(N):
            m = 1
            if j==0 or j==N-1:
                m = sqrt(1.0/2)
            if i==0 or i==N-1:
                m*= sqrt(1./2)
            M[i,j] = m*cos(pi*i*j/(N-1))
    M *= sqrt(2.0 / (N-1))
    return M

def dst1_matrix(N):
    M = zeros((N,N))
    for i in range(N):
        for j in range(N):
            M[i,j] = sin(pi*(i+1)*(j+1)/(N+1))
    M *= sqrt(2.0 / (N+1))
    return M

N=10
M = dct1_matrix(N)
a = arange(N)+1
b = M.dot(a)
c = M.dot(b)
d = dct(a,type=1,norm='ortho')
e = idct(d,type=1,norm='ortho')
print(a)
print(b)
print(c)
print(d)
print(e)
print("dct1 specrta close: ",allclose(b,d))
print("dct1 inverse close: ",allclose(c,e))
print("dct1 original close:",allclose(e,a))

M = dst1_matrix(N)
a = arange(N)+1
b = M.dot(a)
c = M.dot(b)
d = dst(a,type=1,norm='ortho')
e = idst(d,type=1,norm='ortho')
print(a)
print(b)
print(c)
print(d)
print(e)
print("dst1 specrta close: ",allclose(b,d))
print("dst1 inverse close: ",allclose(c,e))
print("dst1 original close:",allclose(e,a))

[zywina]

I have decided to expand this PR to include the DCT-IV and DST-IV algorithms. The DCT-IV variant in particular has a number of important applications and is long overdue to be included in scipy.fftpack. Including them completes the set of "even" trigonometric transforms numbered 1 to 4.

The algorithm used is from this paper. The algorithm relies on the DCT-II with some pre and post operations on the data. It shares the same property as the rest of the DCT/DST variants of supporting all input sizes. On my machine, benchmarks have shown it to be ~1% slower than DCT-II.

The code has been tested for correctness against FFTW in unnormalized mode. With orthonormalization, I have added unit tests to check against naive implementations of the algorithms. In all of the affected cases (DCT/DST-1/4) the transforms serve as their own inverse when norm='ortho' is selected, which in itself is a worthwhile unit test.

[larsoner]

I am a bit out of my depth on the c changes but the tests look reasonable to me.

Can you check how long it takes to run the dct tests before and after this PR? The SciPy test suite is already quite demanding so we want to make sure we minimize added test time.

[larsoner]

This should not be indented, it will not render properly:

https://7334-1460385-gh.circle-artifacts.com/0/html-scipyorg/generated/scipy.fftpack.dct.html#scipy.fftpack.dct

[larsoner]

Is this worth putting in a .. note:: directive?

[larsoner]

It would probably be better to format these with something like:

Type IV
^^^^^^^

to let Sphinx do the sub-heading formatting rather than forcing bold. But this is not your fault (you are just following existing convention) so feel free to change it or not

[zywina]

@larsoner, thanks for detailed comments.

Run time of unit tests:

$ python runtests.py -s fftpack
before:
172 passed, 1 xfailed in 4.57 seconds
after:
211 passed, 1 xfailed in 9.86 seconds

Most of the additional runtime is no-doubt the naive implementation. I could remove the unnormalized version of that to save a few seconds. Since we are already comparing unnormalized versus FFTW, that one can be let go without much worry.

[larsoner]

That sounds worth it for the time savings.

[larsoner]

Changes look good to me and output renders nicely:

https://7360-1460385-gh.circle-artifacts.com/0/html-scipyorg/generated/scipy.fftpack.dst.html#scipy.fftpack.dst

Anyone else want to have a look? @zywina if nobody else has comments in a week or so, feel free to ping me for merge

[larsoner]

(@ilayn it looks like you had comments previously)

[zywina]

@larsoner I think we're good to go.

[larsoner]

It says there is a merge conflict, can you resolve it?

[larsoner]

Thanks @zywina!

[rgommers]

Nice. @zywina could you please add a short description to the release notes at https://github.com/scipy/scipy/wiki/Release-note-entries-for-SciPy-1.2.0?