matrix_fft.c

Go to the documentation of this file.

/*
 * This work is licensed under a Creative Commons 
 * Attribution-Noncommercial-Share Alike 3.0 United States License.
 * 
 *    http://creativecommons.org/licenses/by-nc-sa/3.0/us/
 * 
 * You are free:
 * 
 *    to Share - to copy, distribute, display, and perform the work
 *    to Remix - to make derivative works
 * 
 * Under the following conditions:
 * 
 *    Attribution. You must attribute the work in the manner specified by the
 *    author or licensor (but not in any way that suggests that they endorse you
 *    or your use of the work).
 * 
 *    Noncommercial. You may not use this work for commercial purposes.
 * 
 *    Share Alike. If you alter, transform, or build upon this work, you may
 *    distribute the resulting work only under the same or similar license to
 *    this one.
 * 
 * For any reuse or distribution, you must make clear to others the license
 * terms of this work. The best way to do this is by including this header.
 * 
 * Any of the above conditions can be waived if you get permission from the
 * copyright holder.
 * 
 * Apart from the remix rights granted under this license, nothing in this
 * license impairs or restricts the author's moral rights.
 */


#include <jwsc/config.h>

#include <stdlib.h>
#include <inttypes.h>
#include <assert.h>
#include <math.h>
#ifdef JWSC_HAVE_FFTW3_H
#include <fftw3.h>
#endif

#include "jwsc/matrix/matrix.h"
#include "jwsc/matrix/matrix_math.h"


#ifdef JWSC_HAVE_FFTW3F
void real_to_complex_matrix_dft_f
(
    Matrix_cf**     m_out,
    const Matrix_f* m_in
)
{
    uint32_t    num_rows, num_cols;
    Matrix_f*   m_in_copy;
    fftwf_plan  plan;

    num_rows = m_in->num_rows;
    num_cols = m_in->num_cols;

    create_matrix_cf(m_out, num_rows, (num_cols/2)+1);

    /* For fftw planning, input array does not need to be initialized, hence we
     * allow fftw to destroy the input to (possibly) get better efficiency. */
    m_in_copy = NULL;
    create_matrix_f(&m_in_copy, num_rows, num_cols);

    assert((plan = fftwf_plan_dft_r2c_2d(num_rows, num_cols,
                    *(m_in_copy->elts), (fftwf_complex*)*((*m_out)->elts), 
                    FFTW_MEASURE | FFTW_DESTROY_INPUT)) != NULL);

    /* Copy the input, now that planning is over, and execute. */
    copy_matrix_f(&m_in_copy, m_in);
    fftwf_execute(plan);

    fftwf_destroy_plan(plan);
    free_matrix_f(m_in_copy);
}
#endif

#ifdef JWSC_HAVE_FFTW3
void real_to_complex_matrix_dft_d
(
    Matrix_cd**     m_out,
    const Matrix_d* m_in
)
{
    uint32_t    num_rows, num_cols;
    Matrix_d*   m_in_copy;
    fftw_plan   plan;

    num_rows = m_in->num_rows;
    num_cols = m_in->num_cols;

    create_matrix_cd(m_out, num_rows, (num_cols/2)+1);

    /* For fftw planning, input array does not need to be initialized, hence we
     * allow fftw to destroy the input to (possibly) get better efficiency. */
    m_in_copy = NULL;
    create_matrix_d(&m_in_copy, num_rows, num_cols);

    assert((plan = fftw_plan_dft_r2c_2d(num_rows, num_cols,
                    *(m_in_copy->elts), (fftw_complex*)*((*m_out)->elts), 
                    FFTW_MEASURE | FFTW_DESTROY_INPUT)) != NULL);

    /* Copy the input, now that planning is over, and execute. */
    copy_matrix_d(&m_in_copy, m_in);
    fftw_execute(plan);

    fftw_destroy_plan(plan);
    free_matrix_d(m_in_copy);
}
#endif

#ifdef JWSC_HAVE_FFTW3F
void complex_to_real_matrix_dft_f
(
    Matrix_f**       m_out,
    const Matrix_cf* m_in,
    uint8_t          ncols_even
)
{
    uint32_t   num_rows, num_cols;
    Matrix_cf* m_in_copy;
    fftwf_plan plan;

    num_rows = m_in->num_rows;
    num_cols = (ncols_even) ? (m_in->num_cols-1)*2 : (m_in->num_cols-1)*2 + 1;

    create_matrix_f(m_out, num_rows, num_cols);

    /* For fftw planning, input array does not need to be initialized, hence we
     * allow fftw to destroy the input to (possibly) get better efficiency. */
    m_in_copy = NULL;
    create_matrix_cf(&m_in_copy, m_in->num_rows, m_in->num_cols);

    assert((plan = fftwf_plan_dft_c2r_2d(num_rows, num_cols,
                    (fftwf_complex*)*(m_in_copy->elts), *((*m_out)->elts), 
                    FFTW_MEASURE | FFTW_DESTROY_INPUT)) != NULL);

    /* Copy the input, now that planning is over, and execute. */
    copy_matrix_cf(&m_in_copy, m_in);
    fftwf_execute(plan);

    fftwf_destroy_plan(plan);
    free_matrix_cf(m_in_copy);
}
#endif

#ifdef JWSC_HAVE_FFTW3
void complex_to_real_matrix_dft_d
(
    Matrix_d**       m_out,
    const Matrix_cd* m_in,
    uint8_t          ncols_even
)
{
    uint32_t   num_rows, num_cols;
    Matrix_cd* m_in_copy;
    fftw_plan  plan;

    num_rows = m_in->num_rows;
    num_cols = (ncols_even) ? (m_in->num_cols-1)*2 : (m_in->num_cols-1)*2 + 1;

    create_matrix_d(m_out, num_rows, num_cols);

    /* For fftw planning, input array does not need to be initialized, hence we
     * allow fftw to destroy the input to (possibly) get better efficiency. */
    m_in_copy = NULL;
    create_matrix_cd(&m_in_copy, m_in->num_rows, m_in->num_cols);

    assert((plan = fftw_plan_dft_c2r_2d(num_rows, num_cols,
                    (fftw_complex*)*(m_in_copy->elts), *((*m_out)->elts), 
                    FFTW_MEASURE | FFTW_DESTROY_INPUT)) != NULL);

    /* Copy the input, now that planning is over, and execute. */
    copy_matrix_cd(&m_in_copy, m_in);
    fftw_execute(plan);

    fftw_destroy_plan(plan);
    free_matrix_cd(m_in_copy);
}
#endif

void expand_symmetric_matrix_dft_f
(
    Matrix_cf**      m_out, 
    const Matrix_cf* m_in, 
    uint8_t          ncols_even
)
{
    uint32_t   num_rows, num_cols;
    uint32_t   row, col;
    Matrix_cf* m;

    num_rows = m_in->num_rows;
    num_cols = m_in->num_cols;

    m = (*m_out == m_in) ? NULL : *m_out;
    if (ncols_even)
    {
        create_matrix_cf(&m, num_rows, (m_in->num_cols-1)*2);
    }
    else
    {
        create_matrix_cf(&m, num_rows, (m_in->num_cols-1)*2 + 1);
    }

    assert(copy_matrix_block_into_matrix_cf(m, 0, 0, m_in, 0, 0, 
                m_in->num_rows, m_in->num_cols) == NULL);

    for (row = 0; row < num_rows; row++)
    {
        for (col = m->num_cols-1; col >= m_in->num_cols; col--)
        {
            m->elts[row][col] = m_in->elts[row][m->num_cols-col-1];
        }
    }

    if (*m_out == m_in)
    {
        copy_matrix_cf(m_out, m);
        free_matrix_cf(m);
    }
    else
    {
        *m_out = m;
    }
}

void expand_symmetric_matrix_dft_d
(
    Matrix_cd**      m_out, 
    const Matrix_cd* m_in, 
    uint8_t          ncols_even
)
{
    uint32_t   num_rows, num_cols;
    uint32_t   row, col;
    Matrix_cd* m;

    num_rows = m_in->num_rows;
    num_cols = m_in->num_cols;

    m = (*m_out == m_in) ? NULL : *m_out;
    if (ncols_even)
    {
        create_matrix_cd(&m, num_rows, (m_in->num_cols-1)*2);
    }
    else
    {
        create_matrix_cd(&m, num_rows, (m_in->num_cols-1)*2 + 1);
    }

    assert(copy_matrix_block_into_matrix_cd(m, 0, 0, m_in, 0, 0, 
                m_in->num_rows, m_in->num_cols) == NULL);

    for (row = 0; row < num_rows; row++)
    {
        for (col = m->num_cols-1; col >= m_in->num_cols; col--)
        {
            m->elts[row][col] = m_in->elts[row][m->num_cols-col-1];
        }
    }

    if (*m_out == m_in)
    {
        copy_matrix_cd(m_out, m);
        free_matrix_cd(m);
    }
    else
    {
        *m_out = m;
    }
}

void shrink_symmetric_matrix_dft_f(Matrix_cf** m_out, const Matrix_cf* m_in)
{
    Matrix_cf* m;

    m = (*m_out == m_in) ? NULL : *m_out;
    create_matrix_cf(&m, m_in->num_rows, (m_in->num_cols/2)+1);

    assert(copy_matrix_block_into_matrix_cf(m, 0, 0, m_in, 0, 0, 
                m->num_rows, m->num_cols) == NULL);

    if (*m_out == m_in)
    {
        copy_matrix_cf(m_out, m);
        free_matrix_cf(m);
    }
    else
    {
        *m_out = m;
    }
}

void shrink_symmetric_matrix_dft_d(Matrix_cd** m_out, const Matrix_cd* m_in)
{
    Matrix_cd* m;

    m = (*m_out == m_in) ? NULL : *m_out;
    create_matrix_cd(&m, m_in->num_rows, (m_in->num_cols/2)+1);

    assert(copy_matrix_block_into_matrix_cd(m, 0, 0, m_in, 0, 0, 
                m->num_rows, m->num_cols) == NULL);

    if (*m_out == m_in)
    {
        copy_matrix_cd(m_out, m);
        free_matrix_cd(m);
    }
    else
    {
        *m_out = m;
    }
}

void scale_matrix_dft_f(Matrix_f** m_out, const Matrix_f* m_in)
{
    float scale;

    scale = 1.0f / (float)(m_in->num_rows*m_in->num_cols);
    multiply_matrix_by_scalar_f(m_out, m_in, scale);
}

void scale_matrix_dft_d(Matrix_d** m_out, const Matrix_d* m_in)
{
    float scale;

    scale = 1.0 / (double)(m_in->num_rows*m_in->num_cols);
    multiply_matrix_by_scalar_d(m_out, m_in, scale);
}

void shift_matrix_f
(
    Matrix_f**       m_out,
    const Matrix_f*  m_in
)
{
    uint32_t num_rows, num_cols;
    uint32_t row, col;
    uint32_t half_rows, half_cols;

    Matrix_f* m = NULL;

    num_rows  = m_in->num_rows;
    num_cols  = m_in->num_cols;
    half_rows = num_rows / 2;
    half_cols = num_cols / 2;

    m = (*m_out == m_in) ? NULL : *m_out;
    create_matrix_f(&m, num_rows, num_cols);

    for (row = 0; row < num_rows; row++)
    {
        if (row < half_rows + (num_rows % 2))
        {
            for (col = 0; col < num_cols; col++)
            {
                if (col < half_cols + (num_cols % 2))
                {
                    m->elts[ row ][ col ] = m_in->elts
                        [ row + half_rows ]
                        [ col + half_cols ];
                }
                else
                {
                    m->elts[ row ][ col ] = m_in->elts
                        [ row + half_rows ]
                        [ col-half_cols-(num_cols % 2) ];
                }
            }
        }
        else
        {
            for (col = 0; col < num_cols; col++)
            {
                if (col < half_cols + (num_cols % 2))
                {
                    m->elts[ row ][ col ] = m_in->elts
                        [ row - half_rows - (num_rows % 2) ]
                        [ col + half_cols ];
                }
                else
                {
                    m->elts[ row ][ col ] = m_in->elts
                        [ row - half_rows - (num_rows % 2) ]
                        [ col - half_cols - (num_cols % 2) ];
                }
            }
        }
    }

    if (*m_out == m_in)
    {
        copy_matrix_f(m_out, m);
        free_matrix_f(m);
    }
    else
    {
        *m_out = m;
    }
}

void shift_matrix_d
(
    Matrix_d**       m_out,
    const Matrix_d*  m_in
)
{
    uint32_t num_rows, num_cols;
    uint32_t row, col;
    uint32_t half_rows, half_cols;

    Matrix_d* m = NULL;

    num_rows  = m_in->num_rows;
    num_cols  = m_in->num_cols;
    half_rows = num_rows / 2;
    half_cols = num_cols / 2;

    m = (*m_out == m_in) ? NULL : *m_out;
    create_matrix_d(&m, num_rows, num_cols);

    for (row = 0; row < num_rows; row++)
    {
        if (row < half_rows + (num_rows % 2))
        {
            for (col = 0; col < num_cols; col++)
            {
                if (col < half_cols + (num_cols % 2))
                {
                    m->elts[ row ][ col ] = m_in->elts
                        [ row + half_rows ]
                        [ col + half_cols ];
                }
                else
                {
                    m->elts[ row ][ col ] = m_in->elts
                        [ row + half_rows ]
                        [ col-half_cols-(num_cols % 2) ];
                }
            }
        }
        else
        {
            for (col = 0; col < num_cols; col++)
            {
                if (col < half_cols + (num_cols % 2))
                {
                    m->elts[ row ][ col ] = m_in->elts
                        [ row - half_rows - (num_rows % 2) ]
                        [ col + half_cols ];
                }
                else
                {
                    m->elts[ row ][ col ] = m_in->elts
                        [ row - half_rows - (num_rows % 2) ]
                        [ col - half_cols - (num_cols % 2) ];
                }
            }
        }
    }

    if (*m_out == m_in)
    {
        copy_matrix_d(m_out, m);
        free_matrix_d(m);
    }
    else
    {
        *m_out = m;
    }
}

void shift_matrix_cf
(
    Matrix_cf**       m_out,
    const Matrix_cf*  m_in
)
{
    uint32_t num_rows, num_cols;
    uint32_t row, col;
    uint32_t half_rows, half_cols;

    Matrix_cf* m = NULL;

    num_rows  = m_in->num_rows;
    num_cols  = m_in->num_cols;
    half_rows = num_rows / 2;
    half_cols = num_cols / 2;

    m = (*m_out == m_in) ? NULL : *m_out;
    create_matrix_cf(&m, num_rows, num_cols);

    for (row = 0; row < num_rows; row++)
    {
        if (row < half_rows + (num_rows % 2))
        {
            for (col = 0; col < num_cols; col++)
            {
                if (col < half_cols + (num_cols % 2))
                {
                    m->elts[ row ][ col ] = m_in->elts
                        [ row + half_rows ]
                        [ col + half_cols ];
                }
                else
                {
                    m->elts[ row ][ col ] = m_in->elts
                        [ row + half_rows ]
                        [ col-half_cols-(num_cols % 2) ];
                }
            }
        }
        else
        {
            for (col = 0; col < num_cols; col++)
            {
                if (col < half_cols + (num_cols % 2))
                {
                    m->elts[ row ][ col ] = m_in->elts
                        [ row - half_rows - (num_rows % 2) ]
                        [ col + half_cols ];
                }
                else
                {
                    m->elts[ row ][ col ] = m_in->elts
                        [ row - half_rows - (num_rows % 2) ]
                        [ col - half_cols - (num_cols % 2) ];
                }
            }
        }
    }

    if (*m_out == m_in)
    {
        copy_matrix_cf(m_out, m);
        free_matrix_cf(m);
    }
    else
    {
        *m_out = m;
    }
}

void shift_matrix_cd
(
    Matrix_cd**       m_out,
    const Matrix_cd*  m_in
)
{
    uint32_t num_rows, num_cols;
    uint32_t row, col;
    uint32_t half_rows, half_cols;

    Matrix_cd* m = NULL;

    num_rows  = m_in->num_rows;
    num_cols  = m_in->num_cols;
    half_rows = num_rows / 2;
    half_cols = num_cols / 2;

    m = (*m_out == m_in) ? NULL : *m_out;
    create_matrix_cd(&m, num_rows, num_cols);

    for (row = 0; row < num_rows; row++)
    {
        if (row < half_rows + (num_rows % 2))
        {
            for (col = 0; col < num_cols; col++)
            {
                if (col < half_cols + (num_cols % 2))
                {
                    m->elts[ row ][ col ] = m_in->elts
                        [ row + half_rows ]
                        [ col + half_cols ];
                }
                else
                {
                    m->elts[ row ][ col ] = m_in->elts
                        [ row + half_rows ]
                        [ col-half_cols-(num_cols % 2) ];
                }
            }
        }
        else
        {
            for (col = 0; col < num_cols; col++)
            {
                if (col < half_cols + (num_cols % 2))
                {
                    m->elts[ row ][ col ] = m_in->elts
                        [ row - half_rows - (num_rows % 2) ]
                        [ col + half_cols ];
                }
                else
                {
                    m->elts[ row ][ col ] = m_in->elts
                        [ row - half_rows - (num_rows % 2) ]
                        [ col - half_cols - (num_cols % 2) ];
                }
            }
        }
    }

    if (*m_out == m_in)
    {
        copy_matrix_cd(m_out, m);
        free_matrix_cd(m);
    }
    else
    {
        *m_out = m;
    }
}

void unshift_matrix_f
(
    Matrix_f**       m_out,
    const Matrix_f*  m_in
)
{
    uint32_t num_rows, num_cols;
    uint32_t row, col;
    uint32_t half_rows, half_cols;

    Matrix_f* m = NULL;

    num_rows  = m_in->num_rows;
    num_cols  = m_in->num_cols;
    half_rows = num_rows / 2;
    half_cols = num_cols / 2;

    m = (*m_out == m_in) ? NULL : *m_out;
    create_matrix_f(&m, num_rows, num_cols);

    for (row = 0; row < num_rows; row++)
    {
        if (row < half_rows)
        {
            for (col = 0; col < num_cols; col++)
            {
                if (col < half_cols)
                {
                    m->elts[ row ][ col ] = m_in->elts
                        [ row + half_rows + (num_rows % 2) ]
                        [ col + half_cols + (num_cols % 2) ];
                }
                else
                {
                    m->elts[ row ][ col ] = m_in->elts
                        [ row + half_rows + (num_rows % 2) ]
                        [ col - half_cols ];
                }
            }
        }
        else
        {
            for (col = 0; col < num_cols; col++)
            {
                if (col < half_cols)
                {
                    m->elts[ row ][ col ] = m_in->elts
                        [ row - half_rows ]
                        [ col + half_cols + (num_cols % 2) ];
                }
                else
                {
                    m->elts[ row ][ col ] = m_in->elts
                        [ row - half_rows]
                        [ col - half_cols];
                }
            }
        }
    }

    if (*m_out == m_in)
    {
        copy_matrix_f(m_out, m);
        free_matrix_f(m);
    }
    else
    {
        *m_out = m;
    }
}

void unshift_matrix_d
(
    Matrix_d**       m_out,
    const Matrix_d*  m_in
)
{
    uint32_t num_rows, num_cols;
    uint32_t row, col;
    uint32_t half_rows, half_cols;

    Matrix_d* m = NULL;

    num_rows  = m_in->num_rows;
    num_cols  = m_in->num_cols;
    half_rows = num_rows / 2;
    half_cols = num_cols / 2;

    m = (*m_out == m_in) ? NULL : *m_out;
    create_matrix_d(&m, num_rows, num_cols);

    for (row = 0; row < num_rows; row++)
    {
        if (row < half_rows)
        {
            for (col = 0; col < num_cols; col++)
            {
                if (col < half_cols)
                {
                    m->elts[ row ][ col ] = m_in->elts
                        [ row + half_rows + (num_rows % 2) ]
                        [ col + half_cols + (num_cols % 2) ];
                }
                else
                {
                    m->elts[ row ][ col ] = m_in->elts
                        [ row + half_rows + (num_rows % 2) ]
                        [ col - half_cols ];
                }
            }
        }
        else
        {
            for (col = 0; col < num_cols; col++)
            {
                if (col < half_cols)
                {
                    m->elts[ row ][ col ] = m_in->elts
                        [ row - half_rows ]
                        [ col + half_cols + (num_cols % 2) ];
                }
                else
                {
                    m->elts[ row ][ col ] = m_in->elts
                        [ row - half_rows]
                        [ col - half_cols];
                }
            }
        }
    }

    if (*m_out == m_in)
    {
        copy_matrix_d(m_out, m);
        free_matrix_d(m);
    }
    else
    {
        *m_out = m;
    }
}

void unshift_matrix_cf
(
    Matrix_cf**       m_out,
    const Matrix_cf*  m_in
)
{
    uint32_t num_rows, num_cols;
    uint32_t row, col;
    uint32_t half_rows, half_cols;

    Matrix_cf* m = NULL;

    num_rows  = m_in->num_rows;
    num_cols  = m_in->num_cols;
    half_rows = num_rows / 2;
    half_cols = num_cols / 2;

    m = (*m_out == m_in) ? NULL : *m_out;
    create_matrix_cf(&m, num_rows, num_cols);

    for (row = 0; row < num_rows; row++)
    {
        if (row < half_rows)
        {
            for (col = 0; col < num_cols; col++)
            {
                if (col < half_cols)
                {
                    m->elts[ row ][ col ] = m_in->elts
                        [ row + half_rows + (num_rows % 2) ]
                        [ col + half_cols + (num_cols % 2) ];
                }
                else
                {
                    m->elts[ row ][ col ] = m_in->elts
                        [ row + half_rows + (num_rows % 2) ]
                        [ col - half_cols ];
                }
            }
        }
        else
        {
            for (col = 0; col < num_cols; col++)
            {
                if (col < half_cols)
                {
                    m->elts[ row ][ col ] = m_in->elts
                        [ row - half_rows ]
                        [ col + half_cols + (num_cols % 2) ];
                }
                else
                {
                    m->elts[ row ][ col ] = m_in->elts
                        [ row - half_rows]
                        [ col - half_cols];
                }
            }
        }
    }

    if (*m_out == m_in)
    {
        copy_matrix_cf(m_out, m);
        free_matrix_cf(m);
    }
    else
    {
        *m_out = m;
    }
}

void unshift_matrix_cd
(
    Matrix_cd**       m_out,
    const Matrix_cd*  m_in
)
{
    uint32_t num_rows, num_cols;
    uint32_t row, col;
    uint32_t half_rows, half_cols;

    Matrix_cd* m = NULL;

    num_rows  = m_in->num_rows;
    num_cols  = m_in->num_cols;
    half_rows = num_rows / 2;
    half_cols = num_cols / 2;

    m = (*m_out == m_in) ? NULL : *m_out;
    create_matrix_cd(&m, num_rows, num_cols);

    for (row = 0; row < num_rows; row++)
    {
        if (row < half_rows)
        {
            for (col = 0; col < num_cols; col++)
            {
                if (col < half_cols)
                {
                    m->elts[ row ][ col ] = m_in->elts
                        [ row + half_rows + (num_rows % 2) ]
                        [ col + half_cols + (num_cols % 2) ];
                }
                else
                {
                    m->elts[ row ][ col ] = m_in->elts
                        [ row + half_rows + (num_rows % 2) ]
                        [ col - half_cols ];
                }
            }
        }
        else
        {
            for (col = 0; col < num_cols; col++)
            {
                if (col < half_cols)
                {
                    m->elts[ row ][ col ] = m_in->elts
                        [ row - half_rows ]
                        [ col + half_cols + (num_cols % 2) ];
                }
                else
                {
                    m->elts[ row ][ col ] = m_in->elts
                        [ row - half_rows]
                        [ col - half_cols];
                }
            }
        }
    }

    if (*m_out == m_in)
    {
        copy_matrix_cd(m_out, m);
        free_matrix_cd(m);
    }
    else
    {
        *m_out = m;
    }
}