lbajardsilogic@0: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
lbajardsilogic@0: 
lbajardsilogic@0: /*
lbajardsilogic@0:     Sonic Visualiser
lbajardsilogic@0:     An audio file viewer and annotation editor.
lbajardsilogic@0:     Centre for Digital Music, Queen Mary, University of London.
lbajardsilogic@0:     This file copyright 2006 Chris Cannam and QMUL.
lbajardsilogic@0:     
lbajardsilogic@0:     This program is free software; you can redistribute it and/or
lbajardsilogic@0:     modify it under the terms of the GNU General Public License as
lbajardsilogic@0:     published by the Free Software Foundation; either version 2 of the
lbajardsilogic@0:     License, or (at your option) any later version.  See the file
lbajardsilogic@0:     COPYING included with this distribution for more information.
lbajardsilogic@0: */
lbajardsilogic@0: 
lbajardsilogic@0: #include <iostream>
lbajardsilogic@0: 
lbajardsilogic@0: #include <QMutexLocker>
lbajardsilogic@0: 
lbajardsilogic@191: #include "fileio/MatrixFile.h"
lbajardsilogic@191: #include "base/Profiler.h"
lbajardsilogic@191: 
lbajardsilogic@191: #include "FFTFileCache.h"
lbajardsilogic@191: 
lbajardsilogic@191: 
lbajardsilogic@191: 
lbajardsilogic@191: 
lbajardsilogic@0: // The underlying matrix has height (m_height * 2 + 1).  In each
lbajardsilogic@0: // column we store magnitude at [0], [2] etc and phase at [1], [3]
lbajardsilogic@0: // etc, and then store the normalization factor (maximum magnitude) at
lbajardsilogic@0: // [m_height * 2].
lbajardsilogic@0: 
lbajardsilogic@0: FFTFileCache::FFTFileCache(QString fileBase, MatrixFile::Mode mode,
lbajardsilogic@0:                            StorageType storageType) :
lbajardsilogic@0:     m_writebuf(0),
lbajardsilogic@0:     m_readbuf(0),
lbajardsilogic@0:     m_readbufCol(0),
lbajardsilogic@0:     m_readbufWidth(0),
lbajardsilogic@0:     m_mfc(new MatrixFile
lbajardsilogic@0:           (fileBase, mode, 
lbajardsilogic@0:            storageType == Compact ? sizeof(uint16_t) : sizeof(float),
lbajardsilogic@0:            mode == MatrixFile::ReadOnly)),
lbajardsilogic@0:     m_storageType(storageType)
lbajardsilogic@0: {
lbajardsilogic@0: //    std::cerr << "FFTFileCache: storage type is " << (storageType == Compact ? "Compact" : storageType == Polar ? "Polar" : "Rectangular") << std::endl;
lbajardsilogic@0: }
lbajardsilogic@0: 
lbajardsilogic@0: FFTFileCache::~FFTFileCache()
lbajardsilogic@0: {
lbajardsilogic@0:     if (m_readbuf) delete[] m_readbuf;
lbajardsilogic@0:     if (m_writebuf) delete[] m_writebuf;
lbajardsilogic@0:     delete m_mfc;
lbajardsilogic@0: }
lbajardsilogic@0: 
lbajardsilogic@0: size_t
lbajardsilogic@0: FFTFileCache::getWidth() const
lbajardsilogic@0: {
lbajardsilogic@0:     return m_mfc->getWidth();
lbajardsilogic@0: }
lbajardsilogic@0: 
lbajardsilogic@0: size_t
lbajardsilogic@0: FFTFileCache::getHeight() const
lbajardsilogic@0: {
lbajardsilogic@0:     size_t mh = m_mfc->getHeight();
lbajardsilogic@0:     if (mh > 0) return (mh - 1) / 2;
lbajardsilogic@0:     else return 0;
lbajardsilogic@0: }
lbajardsilogic@0: 
lbajardsilogic@0: void
lbajardsilogic@0: FFTFileCache::resize(size_t width, size_t height)
lbajardsilogic@0: {
lbajardsilogic@0:     QMutexLocker locker(&m_writeMutex);
lbajardsilogic@0: 
lbajardsilogic@0:     m_mfc->resize(width, height * 2 + 1);
lbajardsilogic@0:     if (m_readbuf) {
lbajardsilogic@0:         delete[] m_readbuf;
lbajardsilogic@0:         m_readbuf = 0;
lbajardsilogic@0:     }
lbajardsilogic@0:     if (m_writebuf) {
lbajardsilogic@0:         delete[] m_writebuf;
lbajardsilogic@0:     }
lbajardsilogic@0:     m_writebuf = new char[(height * 2 + 1) * m_mfc->getCellSize()];
lbajardsilogic@0: }
lbajardsilogic@0: 
lbajardsilogic@0: void
lbajardsilogic@0: FFTFileCache::reset()
lbajardsilogic@0: {
lbajardsilogic@0:     m_mfc->reset();
lbajardsilogic@0: }
lbajardsilogic@0: 
lbajardsilogic@0: float
lbajardsilogic@0: FFTFileCache::getMagnitudeAt(size_t x, size_t y) const
lbajardsilogic@0: {
lbajardsilogic@0:     Profiler profiler("FFTFileCache::getMagnitudeAt", false);
lbajardsilogic@0: 
lbajardsilogic@0:     float value = 0.f;
lbajardsilogic@0: 
lbajardsilogic@0:     switch (m_storageType) {
lbajardsilogic@0: 
lbajardsilogic@0:     case Compact:
lbajardsilogic@0:         value = (getFromReadBufCompactUnsigned(x, y * 2) / 65535.0)
lbajardsilogic@0:             * getNormalizationFactor(x);
lbajardsilogic@0:         break;
lbajardsilogic@0: 
lbajardsilogic@0:     case Rectangular:
lbajardsilogic@0:     {
lbajardsilogic@0:         float real, imag;
lbajardsilogic@0:         getValuesAt(x, y, real, imag);
lbajardsilogic@0:         value = sqrtf(real * real + imag * imag);
lbajardsilogic@0:         break;
lbajardsilogic@0:     }
lbajardsilogic@0: 
lbajardsilogic@0:     case Polar:
lbajardsilogic@0:         value = getFromReadBufStandard(x, y * 2);
lbajardsilogic@0:         break;
lbajardsilogic@0:     }
lbajardsilogic@0: 
lbajardsilogic@0:     return value;
lbajardsilogic@0: }
lbajardsilogic@0: 
lbajardsilogic@0: float
lbajardsilogic@0: FFTFileCache::getNormalizedMagnitudeAt(size_t x, size_t y) const
lbajardsilogic@0: {
lbajardsilogic@0:     float value = 0.f;
lbajardsilogic@0: 
lbajardsilogic@0:     switch (m_storageType) {
lbajardsilogic@0: 
lbajardsilogic@0:     case Compact:
lbajardsilogic@0:         value = getFromReadBufCompactUnsigned(x, y * 2) / 65535.0;
lbajardsilogic@0:         break;
lbajardsilogic@0: 
lbajardsilogic@0:     default:
lbajardsilogic@0:     {
lbajardsilogic@0:         float mag = getMagnitudeAt(x, y);
lbajardsilogic@0:         float factor = getNormalizationFactor(x);
lbajardsilogic@0:         if (factor != 0) value = mag / factor;
lbajardsilogic@0:         else value = 0.f;
lbajardsilogic@0:         break;
lbajardsilogic@0:     }
lbajardsilogic@0:     }
lbajardsilogic@0: 
lbajardsilogic@0:     return value;
lbajardsilogic@0: }
lbajardsilogic@0: 
lbajardsilogic@0: float
lbajardsilogic@0: FFTFileCache::getMaximumMagnitudeAt(size_t x) const
lbajardsilogic@0: {
lbajardsilogic@0:     return getNormalizationFactor(x);
lbajardsilogic@0: }
lbajardsilogic@0: 
lbajardsilogic@0: float
lbajardsilogic@0: FFTFileCache::getPhaseAt(size_t x, size_t y) const
lbajardsilogic@0: {
lbajardsilogic@0:     float value = 0.f;
lbajardsilogic@0:     
lbajardsilogic@0:     switch (m_storageType) {
lbajardsilogic@0: 
lbajardsilogic@0:     case Compact:
lbajardsilogic@0:         value = (getFromReadBufCompactSigned(x, y * 2 + 1) / 32767.0) * M_PI;
lbajardsilogic@0:         break;
lbajardsilogic@0: 
lbajardsilogic@0:     case Rectangular:
lbajardsilogic@0:     {
lbajardsilogic@0:         float real, imag;
lbajardsilogic@0:         getValuesAt(x, y, real, imag);
lbajardsilogic@0:         value = princargf(atan2f(imag, real));
lbajardsilogic@0:         break;
lbajardsilogic@0:     }
lbajardsilogic@0: 
lbajardsilogic@0:     case Polar:
lbajardsilogic@0:         value = getFromReadBufStandard(x, y * 2 + 1);
lbajardsilogic@0:         break;
lbajardsilogic@0:     }
lbajardsilogic@0: 
lbajardsilogic@0:     return value;
lbajardsilogic@0: }
lbajardsilogic@0: 
lbajardsilogic@0: void
lbajardsilogic@0: FFTFileCache::getValuesAt(size_t x, size_t y, float &real, float &imag) const
lbajardsilogic@0: {
lbajardsilogic@0:     switch (m_storageType) {
lbajardsilogic@0: 
lbajardsilogic@0:     case Rectangular:
lbajardsilogic@0:         real = getFromReadBufStandard(x, y * 2);
lbajardsilogic@0:         imag = getFromReadBufStandard(x, y * 2 + 1);
lbajardsilogic@0:         return;
lbajardsilogic@0: 
lbajardsilogic@0:     default:
lbajardsilogic@0:         float mag = getMagnitudeAt(x, y);
lbajardsilogic@0:         float phase = getPhaseAt(x, y);
lbajardsilogic@0:         real = mag * cosf(phase);
lbajardsilogic@0:         imag = mag * sinf(phase);
lbajardsilogic@0:         return;
lbajardsilogic@0:     }
lbajardsilogic@0: }
lbajardsilogic@0: 
lbajardsilogic@0: bool
lbajardsilogic@0: FFTFileCache::haveSetColumnAt(size_t x) const
lbajardsilogic@0: {
lbajardsilogic@0:     return m_mfc->haveSetColumnAt(x);
lbajardsilogic@0: }
lbajardsilogic@0: 
lbajardsilogic@0: void
lbajardsilogic@0: FFTFileCache::setColumnAt(size_t x, float *mags, float *phases, float factor)
lbajardsilogic@0: {
lbajardsilogic@0:     QMutexLocker locker(&m_writeMutex);
lbajardsilogic@0: 
lbajardsilogic@0:     size_t h = getHeight();
lbajardsilogic@0: 
lbajardsilogic@0:     switch (m_storageType) {
lbajardsilogic@0: 
lbajardsilogic@0:     case Compact:
lbajardsilogic@0:         for (size_t y = 0; y < h; ++y) {
lbajardsilogic@0:             ((uint16_t *)m_writebuf)[y * 2] = uint16_t((mags[y] / factor) * 65535.0);
lbajardsilogic@0:             ((uint16_t *)m_writebuf)[y * 2 + 1] = uint16_t(int16_t((phases[y] * 32767) / M_PI));
lbajardsilogic@0:         }
lbajardsilogic@0:         break;
lbajardsilogic@0: 
lbajardsilogic@0:     case Rectangular:
lbajardsilogic@0:         for (size_t y = 0; y < h; ++y) {
lbajardsilogic@0:             ((float *)m_writebuf)[y * 2] = mags[y] * cosf(phases[y]);
lbajardsilogic@0:             ((float *)m_writebuf)[y * 2 + 1] = mags[y] * sinf(phases[y]);
lbajardsilogic@0:         }
lbajardsilogic@0:         break;
lbajardsilogic@0: 
lbajardsilogic@0:     case Polar:
lbajardsilogic@0:         for (size_t y = 0; y < h; ++y) {
lbajardsilogic@0:             ((float *)m_writebuf)[y * 2] = mags[y];
lbajardsilogic@0:             ((float *)m_writebuf)[y * 2 + 1] = phases[y];
lbajardsilogic@0:         }
lbajardsilogic@0:         break;
lbajardsilogic@0:     }
lbajardsilogic@0: 
lbajardsilogic@0:     static float maxFactor = 0;
lbajardsilogic@0:     if (factor > maxFactor) maxFactor = factor;
lbajardsilogic@0: //    std::cerr << "Normalization factor: " << factor << ", max " << maxFactor << " (height " << getHeight() << ")" << std::endl;
lbajardsilogic@0: 
lbajardsilogic@0:     if (m_storageType == Compact) {
lbajardsilogic@0:         if (factor < 0.f || factor > 1.f) {
lbajardsilogic@0:             std::cerr << "WARNING: FFTFileCache::setColumnAt: Normalization factor " << factor << " out of range" << std::endl;
lbajardsilogic@0:             if (factor < 0.f) factor = 0.f;
lbajardsilogic@0:             if (factor > 1.f) factor = 1.f;
lbajardsilogic@0:         }
lbajardsilogic@0:         ((uint16_t *)m_writebuf)[h * 2] = (uint16_t)(factor * 65535.0);
lbajardsilogic@0:     } else {
lbajardsilogic@0:         ((float *)m_writebuf)[h * 2] = factor;
lbajardsilogic@0:     }
lbajardsilogic@0:     m_mfc->setColumnAt(x, m_writebuf);
lbajardsilogic@0: }
lbajardsilogic@0: 
lbajardsilogic@0: void
lbajardsilogic@0: FFTFileCache::setColumnAt(size_t x, float *real, float *imag)
lbajardsilogic@0: {
lbajardsilogic@0:     QMutexLocker locker(&m_writeMutex);
lbajardsilogic@0: 
lbajardsilogic@0:     size_t h = getHeight();
lbajardsilogic@0: 
lbajardsilogic@0:     float max = 0.0f;
lbajardsilogic@0: 
lbajardsilogic@0:     switch (m_storageType) {
lbajardsilogic@0: 
lbajardsilogic@0:     case Compact:
lbajardsilogic@0:         for (size_t y = 0; y < h; ++y) {
lbajardsilogic@0:             float mag = sqrtf(real[y] * real[y] + imag[y] * imag[y]);
lbajardsilogic@0:             if (mag > max) max = mag;
lbajardsilogic@0:         }
lbajardsilogic@0:         for (size_t y = 0; y < h; ++y) {
lbajardsilogic@0:             float mag = sqrtf(real[y] * real[y] + imag[y] * imag[y]);
lbajardsilogic@0:             float phase = princargf(atan2f(imag[y], real[y]));
lbajardsilogic@0:             ((uint16_t *)m_writebuf)[y * 2] = uint16_t((mag / max) * 65535.0);
lbajardsilogic@0:             ((uint16_t *)m_writebuf)[y * 2 + 1] = uint16_t(int16_t((phase * 32767) / M_PI));
lbajardsilogic@0:         }
lbajardsilogic@0:         break;
lbajardsilogic@0: 
lbajardsilogic@0:     case Rectangular:
lbajardsilogic@0:         for (size_t y = 0; y < h; ++y) {
lbajardsilogic@0:             ((float *)m_writebuf)[y * 2] = real[y];
lbajardsilogic@0:             ((float *)m_writebuf)[y * 2 + 1] = imag[y];
lbajardsilogic@0:             float mag = sqrtf(real[y] * real[y] + imag[y] * imag[y]);
lbajardsilogic@0:             if (mag > max) max = mag;
lbajardsilogic@0:         }
lbajardsilogic@0:         break;
lbajardsilogic@0: 
lbajardsilogic@0:     case Polar:
lbajardsilogic@0:         for (size_t y = 0; y < h; ++y) {
lbajardsilogic@0:             float mag = sqrtf(real[y] * real[y] + imag[y] * imag[y]);
lbajardsilogic@0:             if (mag > max) max = mag;
lbajardsilogic@0:             ((float *)m_writebuf)[y * 2] = mag;
lbajardsilogic@0:             ((float *)m_writebuf)[y * 2 + 1] = princargf(atan2f(imag[y], real[y]));
lbajardsilogic@0:         }
lbajardsilogic@0:         break;
lbajardsilogic@0:     }
lbajardsilogic@0: 
lbajardsilogic@0:     ((float *)m_writebuf)[h * 2] = max;
lbajardsilogic@0:     m_mfc->setColumnAt(x, m_writebuf);
lbajardsilogic@0: }
lbajardsilogic@0: 
lbajardsilogic@0: size_t
lbajardsilogic@0: FFTFileCache::getCacheSize(size_t width, size_t height, StorageType type)
lbajardsilogic@0: {
lbajardsilogic@0:     return (height * 2 + 1) * width *
lbajardsilogic@0:         (type == Compact ? sizeof(uint16_t) : sizeof(float)) +
lbajardsilogic@0:         2 * sizeof(size_t); // matrix file header size
lbajardsilogic@0: }
lbajardsilogic@0: 
lbajardsilogic@0: void
lbajardsilogic@0: FFTFileCache::populateReadBuf(size_t x) const
lbajardsilogic@0: {
lbajardsilogic@0:     Profiler profiler("FFTFileCache::populateReadBuf", false);
lbajardsilogic@0: 
lbajardsilogic@0:     if (!m_readbuf) {
lbajardsilogic@0:         m_readbuf = new char[m_mfc->getHeight() * 2 * m_mfc->getCellSize()];
lbajardsilogic@0:     }
lbajardsilogic@0:     m_mfc->getColumnAt(x, m_readbuf);
lbajardsilogic@0:     if (m_mfc->haveSetColumnAt(x + 1)) {
lbajardsilogic@0:         m_mfc->getColumnAt
lbajardsilogic@0:             (x + 1, m_readbuf + m_mfc->getCellSize() * m_mfc->getHeight());
lbajardsilogic@0:         m_readbufWidth = 2;
lbajardsilogic@0:     } else {
lbajardsilogic@0:         m_readbufWidth = 1;
lbajardsilogic@0:     }
lbajardsilogic@0:     m_readbufCol = x;
lbajardsilogic@0: }
lbajardsilogic@0: