diff --git a/src/libprojectM/Audio/AudioConstants.hpp b/src/libprojectM/Audio/AudioConstants.hpp index 7e4f18a14c..6a9ff74277 100644 --- a/src/libprojectM/Audio/AudioConstants.hpp +++ b/src/libprojectM/Audio/AudioConstants.hpp @@ -1,10 +1,16 @@ #pragma once +#include + namespace libprojectM { namespace Audio { -static constexpr int WaveformSamples = 576; //!< Number of waveform data samples available for rendering a frame. -static constexpr int SpectrumSamples = 512; //!< Number of spectrum analyzer samples. +static constexpr int AudioBufferSamples = 576; //!< Number of waveform data samples stored in the buffer for analysis. +static constexpr int WaveformSamples = 480; //!< Number of waveform data samples available for rendering a frame. +static constexpr int SpectrumSamples = 512; //!< Number of spectrum analyzer samples. + +using WaveformBuffer = std::array; //!< Buffer with waveform data. Only the first WaveformSamples number of samples are valid. +using SpectrumBuffer = std::array; //!< Buffer with spectrum data. } // namespace Audio } // namespace libprojectM diff --git a/src/libprojectM/Audio/BeatDetect.cpp b/src/libprojectM/Audio/BeatDetect.cpp deleted file mode 100755 index 96eca68dc5..0000000000 --- a/src/libprojectM/Audio/BeatDetect.cpp +++ /dev/null @@ -1,151 +0,0 @@ -/** - * projectM -- Milkdrop-esque visualisation SDK - * Copyright (C)2003-2004 projectM Team - * - * This library is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * This library is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with this library; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - * See 'LICENSE.txt' included within this release - * - */ -/** - * Takes sound data from wherever and returns beat detection values - * Uses statistical Energy-Based methods. Very simple - * - * Some stuff was taken from Frederic Patin's beat-detection article, - * you'll find it online - */ - -#include "BeatDetect.hpp" -#include "PCM.hpp" - -#include -#include - -namespace libprojectM { -namespace Audio { - -BeatDetect::BeatDetect(PCM& _pcm) - : pcm(_pcm) -{ -} - - -auto BeatDetect::Reset() noexcept -> void -{ - this->treb = 0; - this->mid = 0; - this->bass = 0; - this->trebAtt = 0; - this->midAtt = 0; - this->bassAtt = 0; - this->volAtt = 0; - this->volOld = 0; -} - - -auto BeatDetect::GetPCMScale() const noexcept -> float -{ - return beatSensitivity; -} - - -auto BeatDetect::GetFrameAudioData() const -> FrameAudioData -{ - FrameAudioData data{}; - - pcm.GetPcm(data.waveformLeft.data(), CHANNEL_L, WaveformSamples); - pcm.GetPcm(data.waveformRight.data(), CHANNEL_R, WaveformSamples); - pcm.GetSpectrum(data.spectrumLeft.data(), CHANNEL_L, SpectrumSamples); - pcm.GetSpectrum(data.spectrumRight.data(), CHANNEL_R, SpectrumSamples); - - data.vol = vol; - data.volAtt = volAtt; - - data.bass = bass; - data.bassAtt = bassAtt; - data.mid = mid; - data.midAtt = midAtt; - data.treb = treb; - data.trebAtt = trebAtt; - - return data; -} - -auto BeatDetect::CalculateBeatStatistics() -> void -{ - volOld = vol; - - std::array const freqL = - [this]() { - std::array freq{}; - pcm.GetSpectrum(freq.data(), CHANNEL_L, fftLength); - return freq; - }(); - std::array const freqR = - [this]() { - std::array freq{}; - pcm.GetSpectrum(freq.data(), CHANNEL_R, fftLength); - return freq; - }(); - - auto const intensityBetween = [&freqL, &freqR](size_t const from, size_t const to) { - return std::accumulate(&freqL[from], &freqL[to], 0.f) + - std::accumulate(&freqR[from], &freqR[to], 0.f); - }; - - auto const& updateBand = - [](float& band, float& bandAtt, LowPassFilter& bandSmoothed, float const bandIntensity) { - bandSmoothed.Update(bandIntensity); - band = bandIntensity / std::max(0.0001f, bandSmoothed.Get()); - bandAtt = .6f * bandAtt + .4f * band; - bandAtt = std::min(bandAtt, 100.f); - band = std::min(band, 100.f); - }; - - static_assert(fftLength >= 256, "fftLength too small"); - std::array constexpr ranges{0, 3, 23, 255}; - - float const bassIntensity = intensityBetween(ranges[0], ranges[1]); - updateBand(bass, bassAtt, bassSmoothed, bassIntensity); - - float const midIntensity = intensityBetween(ranges[1], ranges[2]); - updateBand(mid, midAtt, midSmoothed, midIntensity); - - float const trebIntensity = intensityBetween(ranges[2], ranges[3]); - updateBand(treb, trebAtt, trebSmoothed, trebIntensity); - - float const volIntensity = bassIntensity + midIntensity + trebIntensity; - updateBand(vol, volAtt, volSmoothed, volIntensity); -} - - - -auto BeatDetect::LowPassFilter::Update(float nextValue) noexcept -> void -{ - m_current -= m_buffer[m_bufferPos] / bufferLength; - m_current += nextValue / bufferLength; - m_buffer[m_bufferPos] = nextValue; - - ++m_bufferPos; - m_bufferPos %= bufferLength; -} - - -auto BeatDetect::LowPassFilter::Get() const noexcept -> float -{ - return m_current; -} - -} // namespace Audio -} // namespace libprojectM diff --git a/src/libprojectM/Audio/BeatDetect.hpp b/src/libprojectM/Audio/BeatDetect.hpp deleted file mode 100755 index c0a8bef612..0000000000 --- a/src/libprojectM/Audio/BeatDetect.hpp +++ /dev/null @@ -1,100 +0,0 @@ -/** - * @file BeatDetect.cpp - * @brief Beat detection class. Takes decompressed sound buffers and returns various characteristics. - * - * projectM -- Milkdrop-esque visualisation SDK - * Copyright (C)2003-2007 projectM Team - * - * This library is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * This library is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with this library; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - * See 'LICENSE.txt' included within this release - * - */ -#pragma once - -#include "FrameAudioData.hpp" -#include "PCM.hpp" - -#include - -namespace libprojectM { -namespace Audio { - -class BeatDetect -{ -public: - // We should probably remove pcm from the constructor, - // and just pass it as an argument to CalculateBeatStatistics. - explicit BeatDetect(PCM& pcm); - - auto Reset() noexcept -> void; - - /** - * Calculates and updates information about the beat - */ - auto CalculateBeatStatistics() -> void; - - // getPCMScale() was added to address https://github.com/projectM-visualizer/projectm/issues/161 - // Returning 1.0 results in using the raw PCM data, which can make the presets look pretty unresponsive - // if the application volume is low. - [[nodiscard]] - auto GetPCMScale() const noexcept -> float; - - /** - * @brief Returns a filled FrameAudioData structure for the current frame. - * @return A FrameAudioData structure with waveform, spectrum beat detection data. - */ - [[nodiscard]] - auto GetFrameAudioData() const -> FrameAudioData; - - float beatSensitivity{1.f}; - - float treb{0.f}; - float mid{0.f}; - float bass{0.f}; - float volOld{0.f}; - - float trebAtt{0.f}; - float midAtt{0.f}; - float bassAtt{0.f}; - float vol{0.f}; - float volAtt{0.f}; - - PCM& pcm; - -private: - class LowPassFilter - { - public: - auto - Update(float nextValue) noexcept -> void; - - [[nodiscard]] auto - Get() const noexcept -> float; - - private: - static size_t constexpr bufferLength{80}; - size_t m_bufferPos{0}; - std::array m_buffer{0.f}; - float m_current{0.f}; - }; - - LowPassFilter bassSmoothed; - LowPassFilter midSmoothed; - LowPassFilter trebSmoothed; - LowPassFilter volSmoothed; -}; - -} // namespace Audio -} // namespace libprojectM diff --git a/src/libprojectM/Audio/CMakeLists.txt b/src/libprojectM/Audio/CMakeLists.txt new file mode 100644 index 0000000000..e4ca0a7136 --- /dev/null +++ b/src/libprojectM/Audio/CMakeLists.txt @@ -0,0 +1,24 @@ + +add_library(Audio OBJECT + AudioConstants.hpp + MilkdropFFT.cpp + MilkdropFFT.hpp + FrameAudioData.cpp + FrameAudioData.hpp + PCM.cpp + PCM.hpp + Loudness.cpp + Loudness.hpp + WaveformAligner.cpp + WaveformAligner.hpp + ) + +target_include_directories(Audio + PUBLIC + ${CMAKE_CURRENT_SOURCE_DIR} + ) + +target_link_libraries(Audio + PUBLIC + libprojectM::API + ) \ No newline at end of file diff --git a/src/libprojectM/Audio/Loudness.cpp b/src/libprojectM/Audio/Loudness.cpp new file mode 100644 index 0000000000..2139c92406 --- /dev/null +++ b/src/libprojectM/Audio/Loudness.cpp @@ -0,0 +1,62 @@ +#include "Loudness.hpp" + +#include + +namespace libprojectM { +namespace Audio { + +Loudness::Loudness(Loudness::Band band) + : m_band(band) +{ +} + +void Loudness::Update(const std::array& spectrumSamples, double secondsSinceLastFrame, uint32_t frame) +{ + SumBand(spectrumSamples); + UpdateBandAverage(secondsSinceLastFrame, frame); +} + +auto Loudness::CurrentRelative() const -> float +{ + return m_currentRelative; +} + +auto Loudness::AverageRelative() const -> float +{ + return m_averageRelative; +} + +void Loudness::SumBand(const std::array& spectrumSamples) +{ + int start = SpectrumSamples * static_cast(m_band) / 6; + int end = SpectrumSamples * (static_cast(m_band) + 1) / 6; + + m_current = 0.0f; + for (int sample = start; sample < end; sample++) + { + m_current += spectrumSamples[sample]; + } +} + +void Loudness::UpdateBandAverage(double secondsSinceLastFrame, uint32_t frame) +{ + float rate = AdjustRateToFps(m_current > m_average ? 0.2f : 0.5f, secondsSinceLastFrame); + m_average = m_average * rate + m_current * (1.0f - rate); + + rate = AdjustRateToFps(frame < 50 ? 0.9f : 0.992f, secondsSinceLastFrame); + m_longAverage = m_longAverage * rate + m_current * (1.0f - rate); + + m_currentRelative = std::fabs(m_longAverage) < 0.001f ? 1.0f : m_current / m_longAverage; + m_averageRelative = std::fabs(m_longAverage) < 0.001f ? 1.0f : m_average / m_longAverage; +} + +auto Loudness::AdjustRateToFps(float rate, double secondsSinceLastFrame) -> float +{ + float const perSecondDecayRateAtFps1 = std::pow(rate, 30.0f); + float const perFrameDecayRateAtFps2 = std::pow(perSecondDecayRateAtFps1, static_cast(secondsSinceLastFrame)); + + return perFrameDecayRateAtFps2; +} + +} // namespace Audio +} // namespace libprojectM diff --git a/src/libprojectM/Audio/Loudness.hpp b/src/libprojectM/Audio/Loudness.hpp new file mode 100644 index 0000000000..55f3baf8cf --- /dev/null +++ b/src/libprojectM/Audio/Loudness.hpp @@ -0,0 +1,97 @@ +/** + * @file Loudness.hpp + * @brief Calculates loudness values in relation to previous frames. +**/ + +#pragma once + +#include "AudioConstants.hpp" + +#include +#include + +namespace libprojectM { +namespace Audio { + +/** + * @brief Calculates beat-detection loudness relative to the previous frame(s). + */ +class Loudness +{ +public: + /** + * @brief Frequency bands. + * Only the first half of the spectrum is used for these bands, each using one third of this half. + */ + enum class Band : int + { + Bass = 0, //!< Bass band (first sixth of the spectrum) + Middles = 1, //!< Middles band (second sixth of the spectrum) + Treble = 2 //!< Treble band (third sixth of the spectrum) + }; + + /** + * @brief Constructor. + * @param band The band to use for this loudness instance. + */ + explicit Loudness(Band band); + + /** + * @brief Updates the beat detection values and averages. + * Must only be called once per frame. + * @param spectrumSamples The current frame's spectrum analyzer samples to use for the update. + * @param secondsSinceLastFrame (Fractional) seconds passed since the last frame. + * @param frame The number of frames already rendered since start. For the first few frames, a different dampening factor + * will be used to avoid "jumpy" behaviour of the values. + */ + void Update(const std::array& spectrumSamples, double secondsSinceLastFrame, uint32_t frame); + + /** + * @brief Returns the current frame's unattenuated loudness relative to the previous frame. + * This value will revolve around 1.0, with <0.7 being very silent and >1.3 very loud audio. + * @return The current frame's unattenuated loudness relative to the previous frame. + */ + auto CurrentRelative() const -> float; + + /** + * @brief Returns the attenuated loudness averaged over the previous frames. + * This value will revolve around 1.0, with <0.7 being very silent and >1.3 very loud audio. + * Does not change as much as the value returned by CurrentRelative(). + * @return The attenuated loudness averaged over the previous frames. + */ + auto AverageRelative() const -> float; + +private: + /** + * @brief Sums up the spectrum samples for the configured band. + * @param spectrumSamples The spectrum analyzer samples to use. + */ + void SumBand(const std::array& spectrumSamples); + + /** + * @brief Updates the short- and long-term averages and relative values. + * @param secondsSinceLastFrame (Fractional) seconds passed since the last frame. + * @param frame The number of frames already rendered since start. + */ + void UpdateBandAverage(double secondsSinceLastFrame, uint32_t frame); + + /** + * @brief Adjusts the dampening rate according the the current FPS. + * @param rate The rate to be dampened. + * @param secondsSinceLastFrame (Fractional) seconds passed since the last frame. + * @return The dampened rate value. + */ + static auto AdjustRateToFps(float rate, double secondsSinceLastFrame) -> float; + + Band m_band{Band::Bass}; //!< The frequency band to use for this instance. + + float m_current{}; //!< The current frame's sum of all frequency strengths in the current band. + float m_average{}; //!< The short-term averaged value of m_current. + float m_longAverage{}; //!< The long-term averaged value of m_current. + + float m_currentRelative{1.0f}; //!< The relative loudness value to the previous frame. + float m_averageRelative{1.0f}; //!< The attenuated relative loudness value. +}; + +} // namespace Audio +} // namespace libprojectM diff --git a/src/libprojectM/Audio/MilkdropFFT.cpp b/src/libprojectM/Audio/MilkdropFFT.cpp new file mode 100644 index 0000000000..83dfe128d4 --- /dev/null +++ b/src/libprojectM/Audio/MilkdropFFT.cpp @@ -0,0 +1,206 @@ +/* + LICENSE + ------- +Copyright 2005-2013 Nullsoft, Inc. +All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, +are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + + * Neither the name of Nullsoft nor the names of its contributors may be used to + endorse or promote products derived from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR +IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR +CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER +IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT +OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +*/ + +#include "MilkdropFFT.hpp" + +namespace libprojectM { +namespace Audio { + +constexpr auto PI = 3.141592653589793238462643383279502884197169399f; + +MilkdropFFT::MilkdropFFT(size_t samplesIn, size_t samplesOut, bool equalize, float envelopePower) + : m_samplesIn(samplesIn) + , m_numFrequencies(samplesOut * 2) +{ + InitBitRevTable(); + InitCosSinTable(); + InitEnvelopeTable(envelopePower); + InitEqualizeTable(equalize); +} + +void MilkdropFFT::InitEnvelopeTable(float power) +{ + if (power < 0.0f) + { + // Keep all values as-is. + m_envelope = std::vector(m_samplesIn, 1.0f); + return; + } + + float const multiplier = 1.0f / static_cast(m_samplesIn) * 2.0f * PI; + + m_envelope.resize(m_samplesIn); + + if (power == 1.0f) + { + for (size_t i = 0; i < m_samplesIn; i++) + { + m_envelope[i] = 0.5f + 0.5f * std::sin(static_cast(i) * multiplier - PI * 0.5f); + } + } + else + { + for (size_t i = 0; i < m_samplesIn; i++) + { + m_envelope[i] = std::pow(0.5f + 0.5f * std::sin(static_cast(i) * multiplier - PI * 0.5f), power); + } + } +} + +void MilkdropFFT::InitEqualizeTable(bool equalize) +{ + if (!equalize) + { + m_equalize = std::vector(m_numFrequencies / 2, 1.0f); + return; + } + + float const scaling = -0.02f; + float const inverseHalfNumFrequencies = 1.0f / static_cast(m_numFrequencies / 2); + + m_equalize.resize(m_numFrequencies / 2); + + for (size_t i = 0; i < m_numFrequencies / 2; i++) + { + m_equalize[i] = scaling * std::log(static_cast(m_numFrequencies / 2 - i) * inverseHalfNumFrequencies); + } +} + +void MilkdropFFT::InitBitRevTable() +{ + m_bitRevTable.resize(m_numFrequencies); + + for (size_t i = 0; i < m_numFrequencies; i++) + { + m_bitRevTable[i] = i; + } + + size_t j{}; + for (size_t i = 0; i < m_numFrequencies; i++) + { + if (j > i) + { + size_t const temp{m_bitRevTable[i]}; + m_bitRevTable[i] = m_bitRevTable[j]; + m_bitRevTable[j] = temp; + } + + size_t m = m_numFrequencies >> 1; + + while (m >= 1 && j >= m) + { + j -= m; + m >>= 1; + } + + j += m; + } +} + +void MilkdropFFT::InitCosSinTable() +{ + size_t tabsize{}; + size_t dftsize{2}; + + while (dftsize <= m_numFrequencies) + { + tabsize++; + dftsize <<= 1; + } + + m_cosSinTable.resize(tabsize); + + dftsize = 2; + size_t index{0}; + while (dftsize <= m_numFrequencies) + { + auto const theta{-2.0f * PI / static_cast(dftsize)}; + m_cosSinTable[index] = std::polar(1.0f, theta); // Radius 1 is the unity circle. + index++; + dftsize <<= 1; + } +} + +void MilkdropFFT::TimeToFrequencyDomain(const std::vector& waveformData, std::vector& spectralData) +{ + if (m_bitRevTable.empty() || m_cosSinTable.empty() || waveformData.size() < m_samplesIn) + { + spectralData.clear(); + return; + } + + // 1. Set up input to the FFT + std::vector> spectrumData(m_numFrequencies, std::complex()); + for (size_t i = 0; i < m_numFrequencies; i++) + { + size_t const idx{m_bitRevTable[i]}; + if (idx < m_samplesIn) + { + spectrumData[i].real(waveformData[idx] * m_envelope[idx]); + } + } + + // 2. Perform FFT + size_t dftSize{2}; + size_t octave{0}; + + while (dftSize <= m_numFrequencies) + { + std::complex w{1.0f, 0.0f}; + std::complex const wp{m_cosSinTable[octave]}; + + size_t const hdftsize{dftSize >> 1}; + + for (size_t m = 0; m < hdftsize; m += 1) + { + for (size_t i = m; i < m_numFrequencies; i += dftSize) + { + size_t const j{i + hdftsize}; + std::complex const tempNum{spectrumData[j] * w}; + spectrumData[j] = spectrumData[i] - tempNum; + spectrumData[i] = spectrumData[i] + tempNum; + } + + w *= wp; + } + + dftSize <<= 1; + octave++; + } + + // 3. Take the magnitude & eventually equalize it (on a log10 scale) for output + spectralData.resize(m_numFrequencies / 2); + for (size_t i = 0; i < m_numFrequencies / 2; i++) + { + spectralData[i] = m_equalize[i] * std::abs(spectrumData[i]); + } +} + +} // namespace Audio +} // namespace libprojectM diff --git a/src/libprojectM/Audio/MilkdropFFT.hpp b/src/libprojectM/Audio/MilkdropFFT.hpp new file mode 100644 index 0000000000..cf2806c798 --- /dev/null +++ b/src/libprojectM/Audio/MilkdropFFT.hpp @@ -0,0 +1,159 @@ +/* + LICENSE + ------- +Copyright 2005-2013 Nullsoft, Inc. +All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, +are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + + * Neither the name of Nullsoft nor the names of its contributors may be used to + endorse or promote products derived from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR +IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR +CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER +IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT +OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +*/ + +#pragma once + +#include +#include + +namespace libprojectM { +namespace Audio { + +/** + * @brief Performs a Fast Fourier Transform on audio sample data. + * Also applies an equalizer pattern with an envelope curve to the resulting data to smooth out + * certain artifacts. + */ +class MilkdropFFT +{ +public: + /** + * Initializes the Fast Fourier Transform. + * @param samplesIn Number of waveform samples which will be fed into the FFT + * @param samplesOut Number of frequency samples generated; MUST BE A POWER OF 2. + * @param equalize true to roughly equalize the magnitude of the basses and trebles, + * false to leave them untouched. + * @param envelopePower Specifies the envelope power. Set to any negative value to disable the envelope. + * See InitEnvelopeTable for more info. + */ + MilkdropFFT(size_t samplesIn, size_t samplesOut, bool equalize = true, float envelopePower = 1.0f); + + /** + * @brief Converts time-domain samples into frequency-domain samples. + * The array lengths are the two parameters to Init(). + * + * The last sample of the output data will represent the frequency + * that is 1/4th of the input sampling rate. For example, + * if the input wave data is sampled at 44,100 Hz, then the last + * sample of the spectral data output will represent the frequency + * 11,025 Hz. The first sample will be 0 Hz; the frequencies of + * the rest of the samples vary linearly in between. + * + * Note that since human hearing is limited to the range 200 - 20,000 + * Hz. 200 is a low bass hum; 20,000 is an ear-piercing high shriek. + * + * Each time the frequency doubles, that sounds like going up an octave. + * That means that the difference between 200 and 300 Hz is FAR more + * than the difference between 5000 and 5100, for example! + * + * So, when trying to analyze bass, you'll want to look at (probably) + * the 200-800 Hz range; whereas for treble, you'll want the 1,400 - + * 11,025 Hz range. + * + * If you want to get 3 bands, try it this way: + * a) 11,025 / 200 = 55.125 + * b) to get the number of octaves between 200 and 11,025 Hz, solve for n: + * 2^n = 55.125 + * n = log 55.125 / log 2 + * n = 5.785 + * c) so each band should represent 5.785/3 = 1.928 octaves; the ranges are: + * 1) 200 - 200*2^1.928 or 200 - 761 Hz + * 2) 200*2^1.928 - 200*2^(1.928*2) or 761 - 2897 Hz + * 3) 200*2^(1.928*2) - 200*2^(1.928*3) or 2897 - 11025 Hz + * + * A simple sine-wave-based envelope is convolved with the waveform + * data before doing the FFT, to emeliorate the bad frequency response + * of a square (i.e. nonexistent) filter. + * + * You might want to slightly damp (blur) the input if your signal isn't + * of a very high quality, to reduce high-frequency noise that would + * otherwise show up in the output. + * @param waveformData The waveform data to convert. Must contain at least the number of elements passed as samplesIn in Init(). + * @param spectralData The resulting frequency data. Vector will be resized to samplesOut elements as passed to Init(). + * If the conversion failed, e.g. not initialized or too few input samples, the result vector will be empty. + */ + void TimeToFrequencyDomain(const std::vector& waveformData, std::vector& spectralData); + + /** + * @brief Returns the number of frequency samples calculated. + * This is twice the value of samplesOut passed to Init(). + * @return The number of frequency samples calculated. + */ + auto NumFrequencies() const -> size_t + { + return m_numFrequencies; + }; + +private: + /** + * @brief Initializes the equalizer envelope table. + * + * This pre-computation is for multiplying the waveform sample + * by a bell-curve-shaped envelope, so we don't see the ugly + * frequency response (oscillations) of a square filter. + * + * - a power of 1.0 will compute the FFT with exactly one convolution. + * - a power of 2.0 is like doing it twice; the resulting frequency + * output will be smoothed out and the peaks will be "fatter". + * - a power of 0.5 is closer to not using an envelope, and you start + * to get the frequency response of the square filter as 'power' + * approaches zero; the peaks get tighter and more precise, but + * you also see small oscillations around their bases. + * @param power Number of convolutions in the envelope. + */ + void InitEnvelopeTable(float power); + + /** + * @brief Calculates equalization factors for each frequency domain. + * + * @param equalize If false, all factors will be set to 1, otherwise will calculate a frequency-based multiplier. + */ + void InitEqualizeTable(bool equalize); + + /** + * @brief Builds the sample lookup table for each octave. + */ + void InitBitRevTable(); + + /** + * @brief Builds a table with the Nth roots of unity inputs for the transform. + */ + void InitCosSinTable(); + + size_t m_samplesIn{}; //!< Number of waveform samples to use for the FFT calculation. + size_t m_numFrequencies{}; //!< Number of frequency samples calculated by the FFT. + + std::vector m_bitRevTable; //!< Index table for frequency-specific waveform data lookups. + std::vector m_envelope; //!< Equalizer envelope table. + std::vector m_equalize; //!< Equalization values. + std::vector> m_cosSinTable; //!< Table with complex polar coordinates for the different frequency domains used in the FFT. +}; + +} // namespace Audio +} // namespace libprojectM diff --git a/src/libprojectM/Audio/PCM.cpp b/src/libprojectM/Audio/PCM.cpp index 8c554d6a32..576e4d3e75 100755 --- a/src/libprojectM/Audio/PCM.cpp +++ b/src/libprojectM/Audio/PCM.cpp @@ -1,260 +1,125 @@ -/** - * @file PCM.cpp - * @brief Takes sound data from wherever and hands it back out. - * - * Returns PCM Data or spectrum data, or the derivative of the PCM data - * - * projectM -- Milkdrop-esque visualisation SDK - * Copyright (C)2003-2004 projectM Team - * - * This library is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * This library is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with this library; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - * See 'LICENSE.txt' included within this release - * - */ - #include "PCM.hpp" -#include "fftsg.h" - -#include -#include #include namespace libprojectM { namespace Audio { -/* - * Here is where we try to do auto volume setting. Doing this here - * means that none of the code downstream (waveforms, beatdetect, etc) needs - * to worry about it. - * - * 1) Don't over react to level changes within a song - * 2) Ignore silence/gaps - * - * I don't know if it's necessary to have both sum and max, but that makes - * it easier to experiment... - */ -auto PCM::AutoLevel::UpdateLevel( - size_t const samples, - double const sum, - double const max) -> double -{ - - // This is an arbitrary number that helps control - // a) how quickly the level can change and - // b) keeps code from being affected by how the caller provides data (lot of short buffers, or fewer long buffers) - size_t constexpr autolevelSegment = 4096; - - if (sum / static_cast(samples) < 0.00001) - { - return m_level; - } - m_levelSum += sum; - m_levelax = std::max(m_levelax, max * 1.02); - m_levelSamples += samples; - if (m_levelSamples >= autolevelSegment || m_l0 <= 0) - { - double const maxRecent = std::max(std::max(m_l0, m_l1), std::max(m_l2, m_levelax)); - m_l0 = m_l1; - m_l1 = m_l2; - m_l2 = m_levelax; - m_levelax *= 0.95; - m_levelSamples = 0; - m_levelSum = 0; - m_level = (m_l0 <= 0) ? maxRecent : m_level * 0.96 + maxRecent * 0.04; - m_level = std::max(m_level, 0.0001); - } - return m_level; -} - - template< - size_t lOffset, - size_t rOffset, - size_t stride, int signalAmplitude, int signalOffset, - class SampleType> -void PCM::AddPcm( + typename SampleType> +void PCM::AddToBuffer( SampleType const* const samples, - size_t const count) + uint32_t channels, + size_t const sampleCount) { - float sum = 0; - float max = 0; - for (size_t i = 0; i < count; i++) + if (channels == 0 || sampleCount == 0) { - size_t const j = (m_start + i) % maxSamples; - m_pcmL[j] = (samples[lOffset + i * stride] - signalOffset) / float(signalAmplitude); - m_pcmR[j] = (samples[rOffset + i * stride] - signalOffset) / float(signalAmplitude); - sum += std::abs(m_pcmL[j]) + std::abs(m_pcmR[j]); - max = std::max(std::max(max, std::abs(m_pcmL[j])), std::abs(m_pcmR[j])); + return; } - m_start = (m_start + count) % maxSamples; - m_level = m_leveler.UpdateLevel(count, sum / 2, max); -} - -void PCM::AddMono(float const* const samples, size_t const count) -{ - AddPcm<0, 0, 1, 1, 0>(samples, count); -} -void PCM::AddMono(uint8_t const* const samples, size_t const count) -{ - AddPcm<0, 0, 1, 128, 0>(samples, count); -} -void PCM::AddMono(int16_t const* const samples, size_t const count) -{ - AddPcm<0, 0, 1, 32768, 0>(samples, count); + for (size_t i = 0; i < sampleCount; i++) + { + size_t const bufferOffset = (m_start + i) % AudioBufferSamples; + m_inputBufferL[bufferOffset] = 128.0f * (static_cast(samples[0 + i * channels]) - float(signalOffset)) / float(signalAmplitude); + if (channels > 1) + { + m_inputBufferR[bufferOffset] = 128.0f * (static_cast(samples[1 + i * channels]) - float(signalOffset)) / float(signalAmplitude); + } + else + { + m_inputBufferR[bufferOffset] = m_inputBufferL[bufferOffset]; + } + } + m_start = (m_start + sampleCount) % AudioBufferSamples; } - -void PCM::AddStereo(float const* const samples, size_t const count) +void PCM::Add(float const* const samples, uint32_t channels, size_t const count) { - AddPcm<0, 1, 2, 1, 0>(samples, count); + AddToBuffer<1, 0>(samples, channels, count); } -void PCM::AddStereo(uint8_t const* const samples, size_t const count) +void PCM::Add(uint8_t const* const samples, uint32_t channels, size_t const count) { - AddPcm<0, 1, 2, 128, 0>(samples, count); + AddToBuffer<128, 128>(samples, channels, count); } -void PCM::AddStereo(int16_t const* const samples, size_t const count) +void PCM::Add(int16_t const* const samples, uint32_t channels, size_t const count) { - AddPcm<0, 1, 2, 32768, 0>(samples, count); + AddToBuffer<0, 32768>(samples, channels, count); } - -template -void Rdft( - int const isgn, - std::array& a, - std::array& ip, - std::array& w) +void PCM::UpdateFrameAudioData(double secondsSinceLastFrame, uint32_t frame) { - // per rdft() documentation - // length of ip >= 2+sqrt(n/2) and length of w == n/2 - // n: length of a, power of 2, n >= 2 - // see fftsg.cpp - static_assert(2 * (ipSize - 2) * (ipSize - 2) >= aSize, - "rdft invariant not preserved: length of ip >= 2+sqrt(n/2)"); - static_assert(2 * wSize == aSize, - "rdft invariant not preserved: length of w == n/2"); - static_assert(aSize >= 2, - "rdft invariant not preserved: n >= 2"); - static_assert((aSize & (aSize - 1)) == 0, - "rdft invariant not preserved: n is power of two"); - rdft(aSize, isgn, a.data(), ip.data(), w.data()); -} + // 1. Copy audio data from input buffer + CopyNewWaveformData(m_inputBufferL, m_waveformL); + CopyNewWaveformData(m_inputBufferR, m_waveformR); + // 2. Update spectrum analyzer data for both channels + UpdateSpectrum(m_waveformL, m_spectrumL); + UpdateSpectrum(m_waveformR, m_spectrumR); -// puts sound data requested at provided pointer -// -// samples is number of PCM samples to return -// returned values are normalized from -1 to 1 + // 3. Align waveforms + m_alignL.Align(m_waveformL); + m_alignR.Align(m_waveformR); -void PCM::GetPcm( - float* const data, - CHANNEL const channel, - size_t const samples) const -{ - assert(channel == 0 || channel == 1); + // 4. Update beat detection values + m_bass.Update(m_spectrumL, secondsSinceLastFrame, frame); + m_middles.Update(m_spectrumL, secondsSinceLastFrame, frame); + m_treble.Update(m_spectrumL, secondsSinceLastFrame, frame); - CopyPcm(data, channel, samples); } - -void PCM::GetSpectrum( - float* const data, - CHANNEL const channel, - size_t const samples) +auto PCM::GetFrameAudioData() const -> FrameAudioData { - assert(channel == 0 || channel == 1); - UpdateFftChannel(channel); + FrameAudioData data{}; - auto const& spectrum = channel == 0 ? m_spectrumL : m_spectrumR; - size_t const count = samples <= fftLength ? samples : fftLength; - for (size_t i = 0; i < count; i++) - { - data[i] = spectrum[i]; - } - for (size_t i = count; i < samples; i++) - { - data[0] = 0; - } -} + std::copy(m_waveformL.begin(), m_waveformL.end(), data.waveformLeft.begin()); + std::copy(m_waveformR.begin(), m_waveformR.end(), data.waveformRight.begin()); + std::copy(m_spectrumL.begin(), m_spectrumL.end(), data.spectrumLeft.begin()); + std::copy(m_spectrumR.begin(), m_spectrumR.end(), data.spectrumRight.begin()); -void PCM::ResetAutoLevel() -{ - m_leveler = AutoLevel(); - m_level = 1.0f; -} + data.bass = m_bass.CurrentRelative(); + data.mid = m_middles.CurrentRelative(); + data.treb = m_treble.CurrentRelative(); -void PCM::UpdateFftChannel(size_t const channel) -{ - assert(channel == 0 || channel == 1); + data.bassAtt = m_bass.AverageRelative(); + data.midAtt = m_middles.AverageRelative(); + data.trebAtt = m_treble.AverageRelative(); - auto& freq = channel == 0 ? m_freqL : m_freqR; - CopyPcm(freq.data(), channel, freq.size()); - Rdft(1, freq, m_ip, m_w); + data.vol = (data.bass + data.mid + data.treb) * 0.333f; + data.volAtt = (data.bassAtt + data.midAtt + data.trebAtt) * 0.333f; - // compute magnitude data (m^2 actually) - auto& spectrum = channel == 0 ? m_spectrumL : m_spectrumR; - for (size_t i = 1; i < fftLength; i++) - { - auto const m2 = static_cast(freq[i * 2] * freq[i * 2] + freq[i * 2 + 1] * freq[i * 2 + 1]); - spectrum[i - 1] = static_cast(i) * m2 / fftLength; - } - spectrum[fftLength - 1] = static_cast(freq[1] * freq[1]); + return data; } -// CPP17: std::clamp -auto Clamp(double const x, double const lo, double const hi) -> double +void PCM::UpdateSpectrum(const WaveformBuffer& waveformData, SpectrumBuffer& spectrumData) { - return x > hi ? hi : x < lo ? lo - : x; -} + std::vector waveformSamples(AudioBufferSamples); + std::vector spectrumValues; -// pull data from circular buffer -void PCM::CopyPcm(float* const to, size_t const channel, size_t const count) const -{ - assert(channel == 0 || channel == 1); - assert(count < maxSamples); - auto const& from = channel == 0 ? m_pcmL : m_pcmR; - const double volume = 1.0 / m_level; - for (size_t i = 0, pos = m_start; i < count; i++) + size_t oldI{0}; + for (size_t i = 0; i < AudioBufferSamples; i++) { - if (pos == 0) - { - pos = maxSamples; - } - to[i] = static_cast(from[--pos] * volume); + // Damp the input into the FFT a bit, to reduce high-frequency noise: + waveformSamples[i] = 0.5f * (waveformData[i] + waveformData[oldI]); + oldI = i; } + + m_fft.TimeToFrequencyDomain(waveformSamples, spectrumValues); + + std::copy(spectrumValues.begin(), spectrumValues.end(), spectrumData.begin()); } -void PCM::CopyPcm(double* const to, size_t const channel, size_t const count) const +void PCM::CopyNewWaveformData(const WaveformBuffer& source, WaveformBuffer& destination) { - assert(channel == 0 || channel == 1); - auto const& from = channel == 0 ? m_pcmL : m_pcmR; - double const volume = 1.0 / m_level; - for (size_t i = 0, pos = m_start; i < count; i++) + auto const bufferStartIndex = m_start.load(); + + for (size_t i = 0; i < AudioBufferSamples; i++) { - if (pos == 0) - { - pos = maxSamples; - } - to[i] = from[--pos] * volume; + destination[i] = source[(bufferStartIndex + i) % AudioBufferSamples]; } } + } // namespace Audio } // namespace libprojectM diff --git a/src/libprojectM/Audio/PCM.hpp b/src/libprojectM/Audio/PCM.hpp index bf46f632b6..c425fc55a8 100755 --- a/src/libprojectM/Audio/PCM.hpp +++ b/src/libprojectM/Audio/PCM.hpp @@ -1,36 +1,21 @@ /** - * projectM -- Milkdrop-esque visualisation SDK - * Copyright (C)2003-2007 projectM Team + * @file PCM.hpp + * @brief Takes sound data from the outside, analyzes it and hands it back out. * - * This library is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * This library is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with this library; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - * See 'LICENSE.txt' included within this release - * - */ -/** - * $Id$ - * - * Encapsulation of raw sound buffer. Used in beat detection - * - * $Log$ - */ + * Returns waveform and spectrum data, as well as relative beat detection values. + **/ #pragma once -#include "projectM-4/projectM_export.h" +#include "AudioConstants.hpp" +#include "FrameAudioData.hpp" +#include "Loudness.hpp" +#include "MilkdropFFT.hpp" +#include "WaveformAligner.hpp" + +#include -#include +#include #include #include @@ -38,124 +23,95 @@ namespace libprojectM { namespace Audio { -// FFT_LENGTH is number of magnitude values available from getSpectrum(). -// Internally this is generated using 2xFFT_LENGTH samples per channel. -size_t constexpr fftLength = 512; - -enum CHANNEL -{ - CHANNEL_L = 0, - CHANNEL_0 = 0, - CHANNEL_R = 1, - CHANNEL_1 = 1 -}; - class PROJECTM_EXPORT PCM { public: - /* maximum number of sound samples that are actually stored. */ - static constexpr size_t maxSamples = 2048; - /** - * @brief Adds a mono pcm buffer to the storage + * @brief Adds new interleaved floating-point PCM data to the buffer. + * Left channel is expected at offset 0, right channel at offset 1. Other channels are ignored. * @param samples The buffer to be added + * @param channels The number of channels in the input data. * @param count The amount of samples in the buffer */ - void AddMono(const float* samples, size_t count); - void AddMono(const uint8_t* samples, size_t count); - void AddMono(const int16_t* samples, size_t count); + void Add(const float* samples, uint32_t channels, size_t count); /** - * @brief Adds a stereo pcm buffer to the storage - * @param samples The buffer to be added. - * The channels are expected to be interleaved, LRLR. - * @param count The amount of samples in each channel (not total samples) + * @brief Adds new mono unsigned 8-bit PCM data to the storage + * Left channel is expected at offset 0, right channel at offset 1. Other channels are ignored. + * @param samples The buffer to be added + * @param channels The number of channels in the input data. + * @param count The amount of samples in the buffer */ - void AddStereo(const float* samples, size_t count); - void AddStereo(const uint8_t* samples, size_t count); - void AddStereo(const int16_t* samples, size_t count); + void Add(const uint8_t* samples, uint32_t channels, size_t count); /** - * PCM data - * The returned data will 'wrap' if more than maxsamples are requested. + * @brief Adds new mono signed 16-bit PCM data to the storage + * Left channel is expected at offset 0, right channel at offset 1. Other channels are ignored. + * @param samples The buffer to be added + * @param channels The number of channels in the input data. + * @param count The amount of samples in the buffer */ - void GetPcm(float* data, CHANNEL channel, size_t samples) const; + void Add(const int16_t* samples, uint32_t channels, size_t count); - /** Spectrum data - * The returned data will be zero padded if more than FFT_LENGTH values are requested + /** + * @brief Updates the internal audio data values for rendering the next frame. + * This method must only be called once per frame, as it does some temporal blending + * and alignments using the previous frame(s) as reference. This function will perform: + * - Passing the waveform data into the spectrum analyzer + * - Aligning waveforms to a best-fit match to the previous frame to produce a calmer waveform shape. + * - Calculating the bass/mid/treb values and their attenuated (time-smoothed) versions. + * + * @param secondsSinceLastFrame Time passed since rendering the last frame. Basically 1.0/FPS. + * @param frame Frames rendered since projectM was started. */ - void GetSpectrum(float* data, CHANNEL channel, size_t samples); + void UpdateFrameAudioData(double secondsSinceLastFrame, uint32_t frame); /** - * @brief Resets the auto leveler state. - * Makes sense after pausing/resuming audio or other interruptions. + * @brief Returns a class holding a copy of the current frame audio data. + * @return A FrameAudioData class with waveform, spectrum and other derived values. */ - void ResetAutoLevel(); + auto GetFrameAudioData() const -> FrameAudioData; -protected: - // CPP20: Could use a struct for the first 5 params to clarify on call site - // together with designated initializers +private: template< - size_t lOffset, - size_t rOffset, - size_t stride, int signalAmplitude, int signalOffset, - class SampleType> - void AddPcm(const SampleType* samples, size_t count); + typename SampleType> + void AddToBuffer(const SampleType* samples, uint32_t channel, size_t sampleCount); + + /** + * Updates FFT data + */ + void UpdateSpectrum(const WaveformBuffer& waveformData, SpectrumBuffer& spectrumData); - // copy data out of the circular PCM buffer - void CopyPcm(float* to, size_t channel, size_t count) const; - void CopyPcm(double* to, size_t channel, size_t count) const; + /** + * Copies data out of the circular input buffer into the per-frame waveform buffer. + */ + void CopyNewWaveformData(const WaveformBuffer& source, WaveformBuffer& destination); - // Updates FFT data - void UpdateFftChannel(size_t channel); + // External input buffer + WaveformBuffer m_inputBufferL{0.f}; //!< Circular buffer for left-channel PCM data. + WaveformBuffer m_inputBufferR{0.f}; //!< Circular buffer for right-channel PCM data. + std::atomic m_start{0}; //!< Circular buffer start index. -private: - // mem-usage: - // pcmd 2x2048*4b = 16K - // vdata 2x512x2*8b = 16K - // spectrum 2x512*4b = 4k - // w = 512*8b = 4k - - // circular PCM buffer - // adjust "volume" of PCM data as we go, this simplifies everything downstream... - // normalize to range [-1.0,1.0] - std::array m_pcmL{0.f}; - std::array m_pcmR{0.f}; - size_t m_start{0}; - size_t m_newSamples{0}; - - // raw FFT data - std::array m_freqL{0.0}; - std::array m_freqR{0.0}; - // magnitude data - std::array m_spectrumL{0.f}; - std::array m_spectrumR{0.f}; - - std::array m_w{0.0}; - std::array m_ip{0}; - - // see https://github.com/projectM-visualizer/projectm/issues/161 - class PROJECTM_EXPORT AutoLevel - { - public: - auto UpdateLevel(size_t samples, double sum, double max) -> double; - - private: - double m_level{0.01}; - // accumulate sample data - size_t m_levelSamples{0}; - double m_levelSum{0.0}; - double m_levelax{0.0}; - double m_l0{-1.0}; - double m_l1{-1.0}; - double m_l2{-1.0}; - }; - - // state for tracking audio level - double m_level{1.f}; - AutoLevel m_leveler{}; + // Frame waveform data + WaveformBuffer m_waveformL{0.f}; //!< Left-channel waveform data, aligned. Only the first WaveformSamples number of samples are valid. + WaveformBuffer m_waveformR{0.f}; //!< Right-channel waveform data, aligned. Only the first WaveformSamples number of samples are valid. + + // Frame spectrum data + SpectrumBuffer m_spectrumL{0.f}; //!< Left-channel spectrum data. + SpectrumBuffer m_spectrumR{0.f}; //!< Right-channel spectrum data. + + MilkdropFFT m_fft{WaveformSamples, SpectrumSamples, true}; //!< Spectrum analyzer instance. + + // Alignment data + WaveformAligner m_alignL; //!< Left-channel waveform alignment. + WaveformAligner m_alignR; //!< Left-channel waveform alignment. + + // Frame beat detection values + Loudness m_bass{Loudness::Band::Bass}; //!< Beat detection/volume for the "bass" band. + Loudness m_middles{Loudness::Band::Middles}; //!< Beat detection/volume for the "middles" band. + Loudness m_treble{Loudness::Band::Treble}; //!< Beat detection/volume for the "treble" band. }; } // namespace Audio diff --git a/src/libprojectM/Audio/WaveformAligner.cpp b/src/libprojectM/Audio/WaveformAligner.cpp new file mode 100644 index 0000000000..737c7890d6 --- /dev/null +++ b/src/libprojectM/Audio/WaveformAligner.cpp @@ -0,0 +1,177 @@ +#include "WaveformAligner.hpp" + +#include +#include + +namespace libprojectM { +namespace Audio { + +WaveformAligner::WaveformAligner() +{ + static const size_t maxOctaves{10}; + static const size_t numOctaves{static_cast(std::floor(std::log(AudioBufferSamples - WaveformSamples) / std::log(2.0f)))}; + m_octaves = numOctaves > maxOctaves ? maxOctaves : numOctaves; + + m_aligmentWeights.resize(m_octaves); + m_firstNonzeroWeights.resize(m_octaves); + m_lastNonzeroWeights.resize(m_octaves); + m_octaveSamples.resize(m_octaves); + m_octaveSampleSpacing.resize(m_octaves); + m_oldWaveformMips.resize(m_octaves); + + m_octaveSamples[0] = AudioBufferSamples; + m_octaveSampleSpacing[0] = AudioBufferSamples - WaveformSamples; + for (size_t octave = 1; octave < m_octaves; octave++) + { + m_octaveSamples[octave] = m_octaveSamples[octave - 1] / 2; + m_octaveSampleSpacing[octave] = m_octaveSampleSpacing[octave - 1] / 2; + } +} + +void WaveformAligner::Align(WaveformBuffer& newWaveform) +{ + if (m_octaves < 4) + { + return; + } + + int alignOffset{}; + + std::vector newWaveformMips(m_octaves, WaveformBuffer()); + std::copy(newWaveform.begin(), newWaveform.end(), newWaveformMips[0].begin()); + + // Calculate mip levels + for (size_t octave = 1; octave < m_octaves; octave++) + { + for (size_t sample = 0; sample < m_octaveSamples[octave]; sample++) + { + newWaveformMips[octave][sample] = 0.5f * (newWaveformMips[octave - 1][sample * 2] + newWaveformMips[octave - 1][sample * 2 + 1]); + } + } + + if (!m_alignWaveReady) + { + m_alignWaveReady = true; + for (size_t octave = 0; octave < m_octaves; octave++) + { + // For example: + // m_octaveSampleSpacing[octave] == 4 + // m_octaveSamples[octave] == 36 + // (so we test 32 samples, w/4 offsets) + size_t const compareSamples = m_octaveSamples[octave] - m_octaveSampleSpacing[octave]; + + for (size_t sample = 0; sample < compareSamples; sample++) + { + auto& tempVal = m_aligmentWeights[octave][sample]; + + // Start with pyramid-shaped PDF, from 0..1..0 + if (sample < compareSamples / 2) + { + tempVal = static_cast(sample * 2) / static_cast(compareSamples); + } + else + { + tempVal = static_cast((compareSamples - 1 - sample) * 2) / static_cast(compareSamples); + } + + // TWEAK how much the center matters, vs. the edges: + tempVal = (tempVal - 0.8f) * 5.0f + 0.8f; + + // Clip + if (tempVal > 1.0f) + { + tempVal = 1.0f; + } + if (tempVal < 0.0f) + { + tempVal = 0.0f; + } + } + + size_t sample{}; + while (m_aligmentWeights[octave][sample] == 0 && sample < compareSamples) + { + sample++; + } + m_firstNonzeroWeights[octave] = sample; + + sample = compareSamples - 1; + while (m_aligmentWeights[octave][sample] == 0 && sample >= 0) + { + sample--; + } + m_lastNonzeroWeights[octave] = sample; + } + } + + int sample1{}; + int sample2{static_cast(m_octaveSampleSpacing[m_octaves - 1])}; + + // Find best match for alignment + for (int octave = static_cast(m_octaves) - 1; octave >= 0; octave--) + { + int lowestErrorOffset{-1}; + float lowestErrorAmount{}; + + for (int sample = sample1; sample < sample2; sample++) + { + float errorSum{}; + + for (size_t i = m_firstNonzeroWeights[octave]; i <= m_lastNonzeroWeights[octave]; i++) + { + errorSum += std::abs((newWaveformMips[octave][i + sample] - m_oldWaveformMips[octave][i + sample]) * m_aligmentWeights[octave][i]); + } + + if (lowestErrorOffset == -1 || errorSum < lowestErrorAmount) + { + lowestErrorOffset = static_cast(sample); + lowestErrorAmount = errorSum; + } + } + + // Now use 'lowestErrorOffset' to guide bounds of search in next octave: + // m_octaveSampleSpacing[octave] == 8 + // m_octaveSamples[octave] == 72 + // -say 'lowestErrorOffset' was 2 + // -that corresponds to samples 4 & 5 of the next octave + // -also, expand about this by 2 samples? YES. + // (so we'd test 64 samples, w/8->4 offsets) + if (octave > 0) + { + sample1 = lowestErrorOffset * 2 - 1; + sample2 = lowestErrorOffset * 2 + 2 + 1; + if (sample1 < 0) + { + sample1 = 0; + } + if (sample2 > static_cast(m_octaveSampleSpacing[octave - 1])) + { + sample2 = static_cast(m_octaveSampleSpacing[octave - 1]); + } + } + else + { + alignOffset = lowestErrorOffset; + } + } + + // Store mip levels for the next frame. + m_oldWaveformMips.clear(); + std::copy(newWaveformMips.begin(), newWaveformMips.end(), std::back_inserter(m_oldWaveformMips)); + + // Finally, apply the results by scooting the aligned samples so that they start at index 0. + if (alignOffset > 0) + { + for (size_t sample = 0; sample < WaveformSamples; sample++) + { + newWaveform[sample] = newWaveform[sample + alignOffset]; + } + + // Set remaining samples to zero. + std::fill_n(newWaveform.begin() + WaveformSamples, AudioBufferSamples - WaveformSamples, 0.0f); + } +} + + +} // namespace Audio +} // namespace libprojectM diff --git a/src/libprojectM/Audio/WaveformAligner.hpp b/src/libprojectM/Audio/WaveformAligner.hpp new file mode 100644 index 0000000000..d8102726ae --- /dev/null +++ b/src/libprojectM/Audio/WaveformAligner.hpp @@ -0,0 +1,51 @@ +/** + * @file WaveformAligner.hpp + * @brief Mip-based waveform alignment algorithm + */ + +#pragma once + +#include "AudioConstants.hpp" + +#include +#include + +namespace libprojectM { +namespace Audio { + +/** + * @class WaveformAligner + * @brief Mip-based waveform alignment algorithm + * + * Calculates the absolute error between the previous and current waveforms over several octaves + * and sample offsets, then shifts the new waveform forward to best align with the previous frame. + * This will keep similar features in-place instead of randomly jumping around on each frame and creates + * for a smoother-looking waveform visualization. + */ +class WaveformAligner +{ +public: + WaveformAligner(); + + /** + * @brief Aligns waveforms to a best-fit match to the previous frame. + * @param[in,out] newWaveform The new waveform to be aligned. + */ + void Align(WaveformBuffer& newWaveform); + +private: + bool m_alignWaveReady{false}; //!< Alignment needs special treatment for the first buffer fill. + + std::vector> m_aligmentWeights; //!< Sample weights per octave. + + size_t m_octaves{}; //!< Number of mip-levels/octaves. + std::vector m_octaveSamples; //!< Samples per octave. + std::vector m_octaveSampleSpacing; //!< Space between samples per octave. + + std::vector m_oldWaveformMips; //!< Mip levels of the previous frame's waveform. + std::vector m_firstNonzeroWeights; //!< First non-zero weight sample index for each octave. + std::vector m_lastNonzeroWeights; //!< Last non-zero weight sample index for each octave. +}; + +} // namespace Audio +} // namespace libprojectM diff --git a/src/libprojectM/Audio/fftsg.cpp b/src/libprojectM/Audio/fftsg.cpp deleted file mode 100755 index 43d7534480..0000000000 --- a/src/libprojectM/Audio/fftsg.cpp +++ /dev/null @@ -1,3314 +0,0 @@ -/* -Fast Fourier/Cosine/Sine Transform - dimension :one - data length :power of 2 - decimation :frequency - radix :split-radix - data :inplace - table :use -functions - cdft: Complex Discrete Fourier Transform - rdft: Real Discrete Fourier Transform - ddct: Discrete Cosine Transform - ddst: Discrete Sine Transform - dfct: Cosine Transform of RDFT (Real Symmetric DFT) - dfst: Sine Transform of RDFT (Real Anti-symmetric DFT) -function prototypes - void cdft(int, int, double *, int *, double *); - void rdft(int, int, double *, int *, double *); - void ddct(int, int, double *, int *, double *); - void ddst(int, int, double *, int *, double *); - void dfct(int, double *, double *, int *, double *); - void dfst(int, double *, double *, int *, double *); -macro definitions - USE_CDFT_PTHREADS : default=not defined - CDFT_THREADS_BEGIN_N : must be >= 512, default=8192 - CDFT_4THREADS_BEGIN_N : must be >= 512, default=65536 - USE_CDFT_WINTHREADS : default=not defined - CDFT_THREADS_BEGIN_N : must be >= 512, default=32768 - CDFT_4THREADS_BEGIN_N : must be >= 512, default=524288 - - --------- Complex DFT (Discrete Fourier Transform) -------- - [definition] - - X[k] = sum_j=0^n-1 x[j]*exp(2*pi*i*j*k/n), 0<=k - X[k] = sum_j=0^n-1 x[j]*exp(-2*pi*i*j*k/n), 0<=k - ip[0] = 0; // first time only - cdft(2*n, 1, a, ip, w); - - ip[0] = 0; // first time only - cdft(2*n, -1, a, ip, w); - [parameters] - 2*n :data length (int) - n >= 1, n = power of 2 - a[0...2*n-1] :input/output data (double *) - input data - a[2*j] = Re(x[j]), - a[2*j+1] = Im(x[j]), 0<=j= 2+sqrt(n) - strictly, - length of ip >= - 2+(1<<(int)(log(n+0.5)/log(2))/2). - ip[0],ip[1] are pointers of the cos/sin table. - w[0...n/2-1] :cos/sin table (double *) - w[],ip[] are initialized if ip[0] == 0. - [remark] - Inverse of - cdft(2*n, -1, a, ip, w); - is - cdft(2*n, 1, a, ip, w); - for (j = 0; j <= 2 * n - 1; j++) { - a[j] *= 1.0 / n; - } - . - - --------- Real DFT / Inverse of Real DFT -------- - [definition] - RDFT - R[k] = sum_j=0^n-1 a[j]*cos(2*pi*j*k/n), 0<=k<=n/2 - I[k] = sum_j=0^n-1 a[j]*sin(2*pi*j*k/n), 0 IRDFT (excluding scale) - a[k] = (R[0] + R[n/2]*cos(pi*k))/2 + - sum_j=1^n/2-1 R[j]*cos(2*pi*j*k/n) + - sum_j=1^n/2-1 I[j]*sin(2*pi*j*k/n), 0<=k - ip[0] = 0; // first time only - rdft(n, 1, a, ip, w); - - ip[0] = 0; // first time only - rdft(n, -1, a, ip, w); - [parameters] - n :data length (int) - n >= 2, n = power of 2 - a[0...n-1] :input/output data (double *) - - output data - a[2*k] = R[k], 0<=k - input data - a[2*j] = R[j], 0<=j= 2+sqrt(n/2) - strictly, - length of ip >= - 2+(1<<(int)(log(n/2+0.5)/log(2))/2). - ip[0],ip[1] are pointers of the cos/sin table. - w[0...n/2-1] :cos/sin table (double *) - w[],ip[] are initialized if ip[0] == 0. - [remark] - Inverse of - rdft(n, 1, a, ip, w); - is - rdft(n, -1, a, ip, w); - for (j = 0; j <= n - 1; j++) { - a[j] *= 2.0 / n; - } - . - - --------- DCT (Discrete Cosine Transform) / Inverse of DCT -------- - [definition] - IDCT (excluding scale) - C[k] = sum_j=0^n-1 a[j]*cos(pi*j*(k+1/2)/n), 0<=k DCT - C[k] = sum_j=0^n-1 a[j]*cos(pi*(j+1/2)*k/n), 0<=k - ip[0] = 0; // first time only - ddct(n, 1, a, ip, w); - - ip[0] = 0; // first time only - ddct(n, -1, a, ip, w); - [parameters] - n :data length (int) - n >= 2, n = power of 2 - a[0...n-1] :input/output data (double *) - output data - a[k] = C[k], 0<=k= 2+sqrt(n/2) - strictly, - length of ip >= - 2+(1<<(int)(log(n/2+0.5)/log(2))/2). - ip[0],ip[1] are pointers of the cos/sin table. - w[0...n*5/4-1] :cos/sin table (double *) - w[],ip[] are initialized if ip[0] == 0. - [remark] - Inverse of - ddct(n, -1, a, ip, w); - is - a[0] *= 0.5; - ddct(n, 1, a, ip, w); - for (j = 0; j <= n - 1; j++) { - a[j] *= 2.0 / n; - } - . - - --------- DST (Discrete Sine Transform) / Inverse of DST -------- - [definition] - IDST (excluding scale) - S[k] = sum_j=1^n A[j]*sin(pi*j*(k+1/2)/n), 0<=k DST - S[k] = sum_j=0^n-1 a[j]*sin(pi*(j+1/2)*k/n), 0 - ip[0] = 0; // first time only - ddst(n, 1, a, ip, w); - - ip[0] = 0; // first time only - ddst(n, -1, a, ip, w); - [parameters] - n :data length (int) - n >= 2, n = power of 2 - a[0...n-1] :input/output data (double *) - - input data - a[j] = A[j], 0 - output data - a[k] = S[k], 0= 2+sqrt(n/2) - strictly, - length of ip >= - 2+(1<<(int)(log(n/2+0.5)/log(2))/2). - ip[0],ip[1] are pointers of the cos/sin table. - w[0...n*5/4-1] :cos/sin table (double *) - w[],ip[] are initialized if ip[0] == 0. - [remark] - Inverse of - ddst(n, -1, a, ip, w); - is - a[0] *= 0.5; - ddst(n, 1, a, ip, w); - for (j = 0; j <= n - 1; j++) { - a[j] *= 2.0 / n; - } - . - - --------- Cosine Transform of RDFT (Real Symmetric DFT) -------- - [definition] - C[k] = sum_j=0^n a[j]*cos(pi*j*k/n), 0<=k<=n - [usage] - ip[0] = 0; // first time only - dfct(n, a, t, ip, w); - [parameters] - n :data length - 1 (int) - n >= 2, n = power of 2 - a[0...n] :input/output data (double *) - output data - a[k] = C[k], 0<=k<=n - t[0...n/2] :work area (double *) - ip[0...*] :work area for bit reversal (int *) - length of ip >= 2+sqrt(n/4) - strictly, - length of ip >= - 2+(1<<(int)(log(n/4+0.5)/log(2))/2). - ip[0],ip[1] are pointers of the cos/sin table. - w[0...n*5/8-1] :cos/sin table (double *) - w[],ip[] are initialized if ip[0] == 0. - [remark] - Inverse of - a[0] *= 0.5; - a[n] *= 0.5; - dfct(n, a, t, ip, w); - is - a[0] *= 0.5; - a[n] *= 0.5; - dfct(n, a, t, ip, w); - for (j = 0; j <= n; j++) { - a[j] *= 2.0 / n; - } - . - - --------- Sine Transform of RDFT (Real Anti-symmetric DFT) -------- - [definition] - S[k] = sum_j=1^n-1 a[j]*sin(pi*j*k/n), 0= 2, n = power of 2 - a[0...n-1] :input/output data (double *) - output data - a[k] = S[k], 0= 2+sqrt(n/4) - strictly, - length of ip >= - 2+(1<<(int)(log(n/4+0.5)/log(2))/2). - ip[0],ip[1] are pointers of the cos/sin table. - w[0...n*5/8-1] :cos/sin table (double *) - w[],ip[] are initialized if ip[0] == 0. - [remark] - Inverse of - dfst(n, a, t, ip, w); - is - dfst(n, a, t, ip, w); - for (j = 1; j <= n - 1; j++) { - a[j] *= 2.0 / n; - } - . - - -Appendix : - The cos/sin table is recalculated when the larger table required. - w[] and ip[] are compatible with all routines. -*/ - - -void cdft(int n, int isgn, double *a, int *ip, double *w) -{ - void makewt(int nw, int *ip, double *w); - void cftfsub(int n, double *a, int *ip, int nw, double *w); - void cftbsub(int n, double *a, int *ip, int nw, double *w); - int nw; - - nw = ip[0]; - if (n > (nw << 2)) { - nw = n >> 2; - makewt(nw, ip, w); - } - if (isgn >= 0) { - cftfsub(n, a, ip, nw, w); - } else { - cftbsub(n, a, ip, nw, w); - } -} - - -void rdft(int n, int isgn, double *a, int *ip, double *w) -{ - void makewt(int nw, int *ip, double *w); - void makect(int nc, int *ip, double *c); - void cftfsub(int n, double *a, int *ip, int nw, double *w); - void cftbsub(int n, double *a, int *ip, int nw, double *w); - void rftfsub(int n, double *a, int nc, double *c); - void rftbsub(int n, double *a, int nc, double *c); - int nw, nc; - double xi; - - nw = ip[0]; - if (n > (nw << 2)) { - nw = n >> 2; - makewt(nw, ip, w); - } - nc = ip[1]; - if (n > (nc << 2)) { - nc = n >> 2; - makect(nc, ip, w + nw); - } - if (isgn >= 0) { - if (n > 4) { - cftfsub(n, a, ip, nw, w); - rftfsub(n, a, nc, w + nw); - } else if (n == 4) { - cftfsub(n, a, ip, nw, w); - } - xi = a[0] - a[1]; - a[0] += a[1]; - a[1] = xi; - } else { - a[1] = 0.5 * (a[0] - a[1]); - a[0] -= a[1]; - if (n > 4) { - rftbsub(n, a, nc, w + nw); - cftbsub(n, a, ip, nw, w); - } else if (n == 4) { - cftbsub(n, a, ip, nw, w); - } - } -} - - -void ddct(int n, int isgn, double *a, int *ip, double *w) -{ - void makewt(int nw, int *ip, double *w); - void makect(int nc, int *ip, double *c); - void cftfsub(int n, double *a, int *ip, int nw, double *w); - void cftbsub(int n, double *a, int *ip, int nw, double *w); - void rftfsub(int n, double *a, int nc, double *c); - void rftbsub(int n, double *a, int nc, double *c); - void dctsub(int n, double *a, int nc, double *c); - int j, nw, nc; - double xr; - - nw = ip[0]; - if (n > (nw << 2)) { - nw = n >> 2; - makewt(nw, ip, w); - } - nc = ip[1]; - if (n > nc) { - nc = n; - makect(nc, ip, w + nw); - } - if (isgn < 0) { - xr = a[n - 1]; - for (j = n - 2; j >= 2; j -= 2) { - a[j + 1] = a[j] - a[j - 1]; - a[j] += a[j - 1]; - } - a[1] = a[0] - xr; - a[0] += xr; - if (n > 4) { - rftbsub(n, a, nc, w + nw); - cftbsub(n, a, ip, nw, w); - } else if (n == 4) { - cftbsub(n, a, ip, nw, w); - } - } - dctsub(n, a, nc, w + nw); - if (isgn >= 0) { - if (n > 4) { - cftfsub(n, a, ip, nw, w); - rftfsub(n, a, nc, w + nw); - } else if (n == 4) { - cftfsub(n, a, ip, nw, w); - } - xr = a[0] - a[1]; - a[0] += a[1]; - for (j = 2; j < n; j += 2) { - a[j - 1] = a[j] - a[j + 1]; - a[j] += a[j + 1]; - } - a[n - 1] = xr; - } -} - - -void ddst(int n, int isgn, double *a, int *ip, double *w) -{ - void makewt(int nw, int *ip, double *w); - void makect(int nc, int *ip, double *c); - void cftfsub(int n, double *a, int *ip, int nw, double *w); - void cftbsub(int n, double *a, int *ip, int nw, double *w); - void rftfsub(int n, double *a, int nc, double *c); - void rftbsub(int n, double *a, int nc, double *c); - void dstsub(int n, double *a, int nc, double *c); - int j, nw, nc; - double xr; - - nw = ip[0]; - if (n > (nw << 2)) { - nw = n >> 2; - makewt(nw, ip, w); - } - nc = ip[1]; - if (n > nc) { - nc = n; - makect(nc, ip, w + nw); - } - if (isgn < 0) { - xr = a[n - 1]; - for (j = n - 2; j >= 2; j -= 2) { - a[j + 1] = -a[j] - a[j - 1]; - a[j] -= a[j - 1]; - } - a[1] = a[0] + xr; - a[0] -= xr; - if (n > 4) { - rftbsub(n, a, nc, w + nw); - cftbsub(n, a, ip, nw, w); - } else if (n == 4) { - cftbsub(n, a, ip, nw, w); - } - } - dstsub(n, a, nc, w + nw); - if (isgn >= 0) { - if (n > 4) { - cftfsub(n, a, ip, nw, w); - rftfsub(n, a, nc, w + nw); - } else if (n == 4) { - cftfsub(n, a, ip, nw, w); - } - xr = a[0] - a[1]; - a[0] += a[1]; - for (j = 2; j < n; j += 2) { - a[j - 1] = -a[j] - a[j + 1]; - a[j] -= a[j + 1]; - } - a[n - 1] = -xr; - } -} - - -void dfct(int n, double *a, double *t, int *ip, double *w) -{ - void makewt(int nw, int *ip, double *w); - void makect(int nc, int *ip, double *c); - void cftfsub(int n, double *a, int *ip, int nw, double *w); - void rftfsub(int n, double *a, int nc, double *c); - void dctsub(int n, double *a, int nc, double *c); - int j, k, l, m, mh, nw, nc; - double xr, xi, yr, yi; - - nw = ip[0]; - if (n > (nw << 3)) { - nw = n >> 3; - makewt(nw, ip, w); - } - nc = ip[1]; - if (n > (nc << 1)) { - nc = n >> 1; - makect(nc, ip, w + nw); - } - m = n >> 1; - yi = a[m]; - xi = a[0] + a[n]; - a[0] -= a[n]; - t[0] = xi - yi; - t[m] = xi + yi; - if (n > 2) { - mh = m >> 1; - for (j = 1; j < mh; j++) { - k = m - j; - xr = a[j] - a[n - j]; - xi = a[j] + a[n - j]; - yr = a[k] - a[n - k]; - yi = a[k] + a[n - k]; - a[j] = xr; - a[k] = yr; - t[j] = xi - yi; - t[k] = xi + yi; - } - t[mh] = a[mh] + a[n - mh]; - a[mh] -= a[n - mh]; - dctsub(m, a, nc, w + nw); - if (m > 4) { - cftfsub(m, a, ip, nw, w); - rftfsub(m, a, nc, w + nw); - } else if (m == 4) { - cftfsub(m, a, ip, nw, w); - } - a[n - 1] = a[0] - a[1]; - a[1] = a[0] + a[1]; - for (j = m - 2; j >= 2; j -= 2) { - a[2 * j + 1] = a[j] + a[j + 1]; - a[2 * j - 1] = a[j] - a[j + 1]; - } - l = 2; - m = mh; - while (m >= 2) { - dctsub(m, t, nc, w + nw); - if (m > 4) { - cftfsub(m, t, ip, nw, w); - rftfsub(m, t, nc, w + nw); - } else if (m == 4) { - cftfsub(m, t, ip, nw, w); - } - a[n - l] = t[0] - t[1]; - a[l] = t[0] + t[1]; - k = 0; - for (j = 2; j < m; j += 2) { - k += l << 2; - a[k - l] = t[j] - t[j + 1]; - a[k + l] = t[j] + t[j + 1]; - } - l <<= 1; - mh = m >> 1; - for (j = 0; j < mh; j++) { - k = m - j; - t[j] = t[m + k] - t[m + j]; - t[k] = t[m + k] + t[m + j]; - } - t[mh] = t[m + mh]; - m = mh; - } - a[l] = t[0]; - a[n] = t[2] - t[1]; - a[0] = t[2] + t[1]; - } else { - a[1] = a[0]; - a[2] = t[0]; - a[0] = t[1]; - } -} - - -void dfst(int n, double *a, double *t, int *ip, double *w) -{ - void makewt(int nw, int *ip, double *w); - void makect(int nc, int *ip, double *c); - void cftfsub(int n, double *a, int *ip, int nw, double *w); - void rftfsub(int n, double *a, int nc, double *c); - void dstsub(int n, double *a, int nc, double *c); - int j, k, l, m, mh, nw, nc; - double xr, xi, yr, yi; - - nw = ip[0]; - if (n > (nw << 3)) { - nw = n >> 3; - makewt(nw, ip, w); - } - nc = ip[1]; - if (n > (nc << 1)) { - nc = n >> 1; - makect(nc, ip, w + nw); - } - if (n > 2) { - m = n >> 1; - mh = m >> 1; - for (j = 1; j < mh; j++) { - k = m - j; - xr = a[j] + a[n - j]; - xi = a[j] - a[n - j]; - yr = a[k] + a[n - k]; - yi = a[k] - a[n - k]; - a[j] = xr; - a[k] = yr; - t[j] = xi + yi; - t[k] = xi - yi; - } - t[0] = a[mh] - a[n - mh]; - a[mh] += a[n - mh]; - a[0] = a[m]; - dstsub(m, a, nc, w + nw); - if (m > 4) { - cftfsub(m, a, ip, nw, w); - rftfsub(m, a, nc, w + nw); - } else if (m == 4) { - cftfsub(m, a, ip, nw, w); - } - a[n - 1] = a[1] - a[0]; - a[1] = a[0] + a[1]; - for (j = m - 2; j >= 2; j -= 2) { - a[2 * j + 1] = a[j] - a[j + 1]; - a[2 * j - 1] = -a[j] - a[j + 1]; - } - l = 2; - m = mh; - while (m >= 2) { - dstsub(m, t, nc, w + nw); - if (m > 4) { - cftfsub(m, t, ip, nw, w); - rftfsub(m, t, nc, w + nw); - } else if (m == 4) { - cftfsub(m, t, ip, nw, w); - } - a[n - l] = t[1] - t[0]; - a[l] = t[0] + t[1]; - k = 0; - for (j = 2; j < m; j += 2) { - k += l << 2; - a[k - l] = -t[j] - t[j + 1]; - a[k + l] = t[j] - t[j + 1]; - } - l <<= 1; - mh = m >> 1; - for (j = 1; j < mh; j++) { - k = m - j; - t[j] = t[m + k] + t[m + j]; - t[k] = t[m + k] - t[m + j]; - } - t[0] = t[m + mh]; - m = mh; - } - a[l] = t[0]; - } - a[0] = 0; -} - - -/* -------- initializing routines -------- */ - - -#include - -void makewt(int nw, int *ip, double *w) -{ - void makeipt(int nw, int *ip); - int j, nwh, nw0, nw1; - double delta, wn4r, wk1r, wk1i, wk3r, wk3i; - - ip[0] = nw; - ip[1] = 1; - if (nw > 2) { - nwh = nw >> 1; - delta = atan(1.0) / nwh; - wn4r = cos(delta * nwh); - w[0] = 1; - w[1] = wn4r; - if (nwh == 4) { - w[2] = cos(delta * 2); - w[3] = sin(delta * 2); - } else if (nwh > 4) { - makeipt(nw, ip); - w[2] = 0.5 / cos(delta * 2); - w[3] = 0.5 / cos(delta * 6); - for (j = 4; j < nwh; j += 4) { - w[j] = cos(delta * j); - w[j + 1] = sin(delta * j); - w[j + 2] = cos(3 * delta * j); - w[j + 3] = -sin(3 * delta * j); - } - } - nw0 = 0; - while (nwh > 2) { - nw1 = nw0 + nwh; - nwh >>= 1; - w[nw1] = 1; - w[nw1 + 1] = wn4r; - if (nwh == 4) { - wk1r = w[nw0 + 4]; - wk1i = w[nw0 + 5]; - w[nw1 + 2] = wk1r; - w[nw1 + 3] = wk1i; - } else if (nwh > 4) { - wk1r = w[nw0 + 4]; - wk3r = w[nw0 + 6]; - w[nw1 + 2] = 0.5 / wk1r; - w[nw1 + 3] = 0.5 / wk3r; - for (j = 4; j < nwh; j += 4) { - wk1r = w[nw0 + 2 * j]; - wk1i = w[nw0 + 2 * j + 1]; - wk3r = w[nw0 + 2 * j + 2]; - wk3i = w[nw0 + 2 * j + 3]; - w[nw1 + j] = wk1r; - w[nw1 + j + 1] = wk1i; - w[nw1 + j + 2] = wk3r; - w[nw1 + j + 3] = wk3i; - } - } - nw0 = nw1; - } - } -} - - -void makeipt(int nw, int *ip) -{ - int j, l, m, m2, p, q; - - ip[2] = 0; - ip[3] = 16; - m = 2; - for (l = nw; l > 32; l >>= 2) { - m2 = m << 1; - q = m2 << 3; - for (j = m; j < m2; j++) { - p = ip[j] << 2; - ip[m + j] = p; - ip[m2 + j] = p + q; - } - m = m2; - } -} - - -void makect(int nc, int *ip, double *c) -{ - int j, nch; - double delta; - - ip[1] = nc; - if (nc > 1) { - nch = nc >> 1; - delta = atan(1.0) / nch; - c[0] = cos(delta * nch); - c[nch] = 0.5 * c[0]; - for (j = 1; j < nch; j++) { - c[j] = 0.5 * cos(delta * j); - c[nc - j] = 0.5 * sin(delta * j); - } - } -} - - -/* -------- child routines -------- */ - - -#ifdef USE_CDFT_PTHREADS -#define USE_CDFT_THREADS -#ifndef CDFT_THREADS_BEGIN_N -#define CDFT_THREADS_BEGIN_N 8192 -#endif -#ifndef CDFT_4THREADS_BEGIN_N -#define CDFT_4THREADS_BEGIN_N 65536 -#endif -#include -#include -#include -#define cdft_thread_t pthread_t -#define cdft_thread_create(thp,func,argp) { \ - if (pthread_create(thp, NULL, func, (void *) argp) != 0) { \ - fprintf(stderr, "cdft thread error\n"); \ - exit(1); \ - } \ -} -#define cdft_thread_wait(th) { \ - if (pthread_join(th, NULL) != 0) { \ - fprintf(stderr, "cdft thread error\n"); \ - exit(1); \ - } \ -} -#endif /* USE_CDFT_PTHREADS */ - - -#ifdef USE_CDFT_WINTHREADS -#define USE_CDFT_THREADS -#ifndef CDFT_THREADS_BEGIN_N -#define CDFT_THREADS_BEGIN_N 32768 -#endif -#ifndef CDFT_4THREADS_BEGIN_N -#define CDFT_4THREADS_BEGIN_N 524288 -#endif -#include -#include -#include -#define cdft_thread_t HANDLE -#define cdft_thread_create(thp,func,argp) { \ - DWORD thid; \ - *(thp) = CreateThread(NULL, 0, (LPTHREAD_START_ROUTINE) func, (LPVOID) argp, 0, &thid); \ - if (*(thp) == 0) { \ - fprintf(stderr, "cdft thread error\n"); \ - exit(1); \ - } \ -} -#define cdft_thread_wait(th) { \ - WaitForSingleObject(th, INFINITE); \ - CloseHandle(th); \ -} -#endif /* USE_CDFT_WINTHREADS */ - - -void cftfsub(int n, double *a, int *ip, int nw, double *w) -{ - void bitrv2(int n, int *ip, double *a); - void bitrv216(double *a); - void bitrv208(double *a); - void cftf1st(int n, double *a, double *w); - void cftrec4(int n, double *a, int nw, double *w); - void cftleaf(int n, int isplt, double *a, int nw, double *w); - void cftfx41(int n, double *a, int nw, double *w); - void cftf161(double *a, double *w); - void cftf081(double *a, double *w); - void cftf040(double *a); - void cftx020(double *a); -#ifdef USE_CDFT_THREADS - void cftrec4_th(int n, double *a, int nw, double *w); -#endif /* USE_CDFT_THREADS */ - - if (n > 8) { - if (n > 32) { - cftf1st(n, a, &w[nw - (n >> 2)]); -#ifdef USE_CDFT_THREADS - if (n > CDFT_THREADS_BEGIN_N) { - cftrec4_th(n, a, nw, w); - } else -#endif /* USE_CDFT_THREADS */ - if (n > 512) { - cftrec4(n, a, nw, w); - } else if (n > 128) { - cftleaf(n, 1, a, nw, w); - } else { - cftfx41(n, a, nw, w); - } - bitrv2(n, ip, a); - } else if (n == 32) { - cftf161(a, &w[nw - 8]); - bitrv216(a); - } else { - cftf081(a, w); - bitrv208(a); - } - } else if (n == 8) { - cftf040(a); - } else if (n == 4) { - cftx020(a); - } -} - - -void cftbsub(int n, double *a, int *ip, int nw, double *w) -{ - void bitrv2conj(int n, int *ip, double *a); - void bitrv216neg(double *a); - void bitrv208neg(double *a); - void cftb1st(int n, double *a, double *w); - void cftrec4(int n, double *a, int nw, double *w); - void cftleaf(int n, int isplt, double *a, int nw, double *w); - void cftfx41(int n, double *a, int nw, double *w); - void cftf161(double *a, double *w); - void cftf081(double *a, double *w); - void cftb040(double *a); - void cftx020(double *a); -#ifdef USE_CDFT_THREADS - void cftrec4_th(int n, double *a, int nw, double *w); -#endif /* USE_CDFT_THREADS */ - - if (n > 8) { - if (n > 32) { - cftb1st(n, a, &w[nw - (n >> 2)]); -#ifdef USE_CDFT_THREADS - if (n > CDFT_THREADS_BEGIN_N) { - cftrec4_th(n, a, nw, w); - } else -#endif /* USE_CDFT_THREADS */ - if (n > 512) { - cftrec4(n, a, nw, w); - } else if (n > 128) { - cftleaf(n, 1, a, nw, w); - } else { - cftfx41(n, a, nw, w); - } - bitrv2conj(n, ip, a); - } else if (n == 32) { - cftf161(a, &w[nw - 8]); - bitrv216neg(a); - } else { - cftf081(a, w); - bitrv208neg(a); - } - } else if (n == 8) { - cftb040(a); - } else if (n == 4) { - cftx020(a); - } -} - - -void bitrv2(int n, int *ip, double *a) -{ - int j, j1, k, k1, l, m, nh, nm; - double xr, xi, yr, yi; - - m = 1; - for (l = n >> 2; l > 8; l >>= 2) { - m <<= 1; - } - nh = n >> 1; - nm = 4 * m; - if (l == 8) { - for (k = 0; k < m; k++) { - for (j = 0; j < k; j++) { - j1 = 4 * j + 2 * ip[m + k]; - k1 = 4 * k + 2 * ip[m + j]; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 += nm; - k1 += 2 * nm; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 += nm; - k1 -= nm; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 += nm; - k1 += 2 * nm; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 += nh; - k1 += 2; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 -= nm; - k1 -= 2 * nm; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 -= nm; - k1 += nm; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 -= nm; - k1 -= 2 * nm; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 += 2; - k1 += nh; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 += nm; - k1 += 2 * nm; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 += nm; - k1 -= nm; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 += nm; - k1 += 2 * nm; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 -= nh; - k1 -= 2; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 -= nm; - k1 -= 2 * nm; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 -= nm; - k1 += nm; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 -= nm; - k1 -= 2 * nm; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - } - k1 = 4 * k + 2 * ip[m + k]; - j1 = k1 + 2; - k1 += nh; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 += nm; - k1 += 2 * nm; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 += nm; - k1 -= nm; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 -= 2; - k1 -= nh; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 += nh + 2; - k1 += nh + 2; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 -= nh - nm; - k1 += 2 * nm - 2; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - } - } else { - for (k = 0; k < m; k++) { - for (j = 0; j < k; j++) { - j1 = 4 * j + ip[m + k]; - k1 = 4 * k + ip[m + j]; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 += nm; - k1 += nm; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 += nh; - k1 += 2; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 -= nm; - k1 -= nm; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 += 2; - k1 += nh; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 += nm; - k1 += nm; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 -= nh; - k1 -= 2; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 -= nm; - k1 -= nm; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - } - k1 = 4 * k + ip[m + k]; - j1 = k1 + 2; - k1 += nh; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 += nm; - k1 += nm; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - } - } -} - - -void bitrv2conj(int n, int *ip, double *a) -{ - int j, j1, k, k1, l, m, nh, nm; - double xr, xi, yr, yi; - - m = 1; - for (l = n >> 2; l > 8; l >>= 2) { - m <<= 1; - } - nh = n >> 1; - nm = 4 * m; - if (l == 8) { - for (k = 0; k < m; k++) { - for (j = 0; j < k; j++) { - j1 = 4 * j + 2 * ip[m + k]; - k1 = 4 * k + 2 * ip[m + j]; - xr = a[j1]; - xi = -a[j1 + 1]; - yr = a[k1]; - yi = -a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 += nm; - k1 += 2 * nm; - xr = a[j1]; - xi = -a[j1 + 1]; - yr = a[k1]; - yi = -a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 += nm; - k1 -= nm; - xr = a[j1]; - xi = -a[j1 + 1]; - yr = a[k1]; - yi = -a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 += nm; - k1 += 2 * nm; - xr = a[j1]; - xi = -a[j1 + 1]; - yr = a[k1]; - yi = -a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 += nh; - k1 += 2; - xr = a[j1]; - xi = -a[j1 + 1]; - yr = a[k1]; - yi = -a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 -= nm; - k1 -= 2 * nm; - xr = a[j1]; - xi = -a[j1 + 1]; - yr = a[k1]; - yi = -a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 -= nm; - k1 += nm; - xr = a[j1]; - xi = -a[j1 + 1]; - yr = a[k1]; - yi = -a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 -= nm; - k1 -= 2 * nm; - xr = a[j1]; - xi = -a[j1 + 1]; - yr = a[k1]; - yi = -a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 += 2; - k1 += nh; - xr = a[j1]; - xi = -a[j1 + 1]; - yr = a[k1]; - yi = -a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 += nm; - k1 += 2 * nm; - xr = a[j1]; - xi = -a[j1 + 1]; - yr = a[k1]; - yi = -a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 += nm; - k1 -= nm; - xr = a[j1]; - xi = -a[j1 + 1]; - yr = a[k1]; - yi = -a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 += nm; - k1 += 2 * nm; - xr = a[j1]; - xi = -a[j1 + 1]; - yr = a[k1]; - yi = -a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 -= nh; - k1 -= 2; - xr = a[j1]; - xi = -a[j1 + 1]; - yr = a[k1]; - yi = -a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 -= nm; - k1 -= 2 * nm; - xr = a[j1]; - xi = -a[j1 + 1]; - yr = a[k1]; - yi = -a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 -= nm; - k1 += nm; - xr = a[j1]; - xi = -a[j1 + 1]; - yr = a[k1]; - yi = -a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 -= nm; - k1 -= 2 * nm; - xr = a[j1]; - xi = -a[j1 + 1]; - yr = a[k1]; - yi = -a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - } - k1 = 4 * k + 2 * ip[m + k]; - j1 = k1 + 2; - k1 += nh; - a[j1 - 1] = -a[j1 - 1]; - xr = a[j1]; - xi = -a[j1 + 1]; - yr = a[k1]; - yi = -a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - a[k1 + 3] = -a[k1 + 3]; - j1 += nm; - k1 += 2 * nm; - xr = a[j1]; - xi = -a[j1 + 1]; - yr = a[k1]; - yi = -a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 += nm; - k1 -= nm; - xr = a[j1]; - xi = -a[j1 + 1]; - yr = a[k1]; - yi = -a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 -= 2; - k1 -= nh; - xr = a[j1]; - xi = -a[j1 + 1]; - yr = a[k1]; - yi = -a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 += nh + 2; - k1 += nh + 2; - xr = a[j1]; - xi = -a[j1 + 1]; - yr = a[k1]; - yi = -a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 -= nh - nm; - k1 += 2 * nm - 2; - a[j1 - 1] = -a[j1 - 1]; - xr = a[j1]; - xi = -a[j1 + 1]; - yr = a[k1]; - yi = -a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - a[k1 + 3] = -a[k1 + 3]; - } - } else { - for (k = 0; k < m; k++) { - for (j = 0; j < k; j++) { - j1 = 4 * j + ip[m + k]; - k1 = 4 * k + ip[m + j]; - xr = a[j1]; - xi = -a[j1 + 1]; - yr = a[k1]; - yi = -a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 += nm; - k1 += nm; - xr = a[j1]; - xi = -a[j1 + 1]; - yr = a[k1]; - yi = -a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 += nh; - k1 += 2; - xr = a[j1]; - xi = -a[j1 + 1]; - yr = a[k1]; - yi = -a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 -= nm; - k1 -= nm; - xr = a[j1]; - xi = -a[j1 + 1]; - yr = a[k1]; - yi = -a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 += 2; - k1 += nh; - xr = a[j1]; - xi = -a[j1 + 1]; - yr = a[k1]; - yi = -a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 += nm; - k1 += nm; - xr = a[j1]; - xi = -a[j1 + 1]; - yr = a[k1]; - yi = -a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 -= nh; - k1 -= 2; - xr = a[j1]; - xi = -a[j1 + 1]; - yr = a[k1]; - yi = -a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 -= nm; - k1 -= nm; - xr = a[j1]; - xi = -a[j1 + 1]; - yr = a[k1]; - yi = -a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - } - k1 = 4 * k + ip[m + k]; - j1 = k1 + 2; - k1 += nh; - a[j1 - 1] = -a[j1 - 1]; - xr = a[j1]; - xi = -a[j1 + 1]; - yr = a[k1]; - yi = -a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - a[k1 + 3] = -a[k1 + 3]; - j1 += nm; - k1 += nm; - a[j1 - 1] = -a[j1 - 1]; - xr = a[j1]; - xi = -a[j1 + 1]; - yr = a[k1]; - yi = -a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - a[k1 + 3] = -a[k1 + 3]; - } - } -} - - -void bitrv216(double *a) -{ - double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i, - x5r, x5i, x7r, x7i, x8r, x8i, x10r, x10i, - x11r, x11i, x12r, x12i, x13r, x13i, x14r, x14i; - - x1r = a[2]; - x1i = a[3]; - x2r = a[4]; - x2i = a[5]; - x3r = a[6]; - x3i = a[7]; - x4r = a[8]; - x4i = a[9]; - x5r = a[10]; - x5i = a[11]; - x7r = a[14]; - x7i = a[15]; - x8r = a[16]; - x8i = a[17]; - x10r = a[20]; - x10i = a[21]; - x11r = a[22]; - x11i = a[23]; - x12r = a[24]; - x12i = a[25]; - x13r = a[26]; - x13i = a[27]; - x14r = a[28]; - x14i = a[29]; - a[2] = x8r; - a[3] = x8i; - a[4] = x4r; - a[5] = x4i; - a[6] = x12r; - a[7] = x12i; - a[8] = x2r; - a[9] = x2i; - a[10] = x10r; - a[11] = x10i; - a[14] = x14r; - a[15] = x14i; - a[16] = x1r; - a[17] = x1i; - a[20] = x5r; - a[21] = x5i; - a[22] = x13r; - a[23] = x13i; - a[24] = x3r; - a[25] = x3i; - a[26] = x11r; - a[27] = x11i; - a[28] = x7r; - a[29] = x7i; -} - - -void bitrv216neg(double *a) -{ - double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i, - x5r, x5i, x6r, x6i, x7r, x7i, x8r, x8i, - x9r, x9i, x10r, x10i, x11r, x11i, x12r, x12i, - x13r, x13i, x14r, x14i, x15r, x15i; - - x1r = a[2]; - x1i = a[3]; - x2r = a[4]; - x2i = a[5]; - x3r = a[6]; - x3i = a[7]; - x4r = a[8]; - x4i = a[9]; - x5r = a[10]; - x5i = a[11]; - x6r = a[12]; - x6i = a[13]; - x7r = a[14]; - x7i = a[15]; - x8r = a[16]; - x8i = a[17]; - x9r = a[18]; - x9i = a[19]; - x10r = a[20]; - x10i = a[21]; - x11r = a[22]; - x11i = a[23]; - x12r = a[24]; - x12i = a[25]; - x13r = a[26]; - x13i = a[27]; - x14r = a[28]; - x14i = a[29]; - x15r = a[30]; - x15i = a[31]; - a[2] = x15r; - a[3] = x15i; - a[4] = x7r; - a[5] = x7i; - a[6] = x11r; - a[7] = x11i; - a[8] = x3r; - a[9] = x3i; - a[10] = x13r; - a[11] = x13i; - a[12] = x5r; - a[13] = x5i; - a[14] = x9r; - a[15] = x9i; - a[16] = x1r; - a[17] = x1i; - a[18] = x14r; - a[19] = x14i; - a[20] = x6r; - a[21] = x6i; - a[22] = x10r; - a[23] = x10i; - a[24] = x2r; - a[25] = x2i; - a[26] = x12r; - a[27] = x12i; - a[28] = x4r; - a[29] = x4i; - a[30] = x8r; - a[31] = x8i; -} - - -void bitrv208(double *a) -{ - double x1r, x1i, x3r, x3i, x4r, x4i, x6r, x6i; - - x1r = a[2]; - x1i = a[3]; - x3r = a[6]; - x3i = a[7]; - x4r = a[8]; - x4i = a[9]; - x6r = a[12]; - x6i = a[13]; - a[2] = x4r; - a[3] = x4i; - a[6] = x6r; - a[7] = x6i; - a[8] = x1r; - a[9] = x1i; - a[12] = x3r; - a[13] = x3i; -} - - -void bitrv208neg(double *a) -{ - double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i, - x5r, x5i, x6r, x6i, x7r, x7i; - - x1r = a[2]; - x1i = a[3]; - x2r = a[4]; - x2i = a[5]; - x3r = a[6]; - x3i = a[7]; - x4r = a[8]; - x4i = a[9]; - x5r = a[10]; - x5i = a[11]; - x6r = a[12]; - x6i = a[13]; - x7r = a[14]; - x7i = a[15]; - a[2] = x7r; - a[3] = x7i; - a[4] = x3r; - a[5] = x3i; - a[6] = x5r; - a[7] = x5i; - a[8] = x1r; - a[9] = x1i; - a[10] = x6r; - a[11] = x6i; - a[12] = x2r; - a[13] = x2i; - a[14] = x4r; - a[15] = x4i; -} - - -void cftf1st(int n, double *a, double *w) -{ - int j, j0, j1, j2, j3, k, m, mh; - double wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i, - wd1r, wd1i, wd3r, wd3i; - double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, - y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i; - - mh = n >> 3; - m = 2 * mh; - j1 = m; - j2 = j1 + m; - j3 = j2 + m; - x0r = a[0] + a[j2]; - x0i = a[1] + a[j2 + 1]; - x1r = a[0] - a[j2]; - x1i = a[1] - a[j2 + 1]; - x2r = a[j1] + a[j3]; - x2i = a[j1 + 1] + a[j3 + 1]; - x3r = a[j1] - a[j3]; - x3i = a[j1 + 1] - a[j3 + 1]; - a[0] = x0r + x2r; - a[1] = x0i + x2i; - a[j1] = x0r - x2r; - a[j1 + 1] = x0i - x2i; - a[j2] = x1r - x3i; - a[j2 + 1] = x1i + x3r; - a[j3] = x1r + x3i; - a[j3 + 1] = x1i - x3r; - wn4r = w[1]; - csc1 = w[2]; - csc3 = w[3]; - wd1r = 1; - wd1i = 0; - wd3r = 1; - wd3i = 0; - k = 0; - for (j = 2; j < mh - 2; j += 4) { - k += 4; - wk1r = csc1 * (wd1r + w[k]); - wk1i = csc1 * (wd1i + w[k + 1]); - wk3r = csc3 * (wd3r + w[k + 2]); - wk3i = csc3 * (wd3i + w[k + 3]); - wd1r = w[k]; - wd1i = w[k + 1]; - wd3r = w[k + 2]; - wd3i = w[k + 3]; - j1 = j + m; - j2 = j1 + m; - j3 = j2 + m; - x0r = a[j] + a[j2]; - x0i = a[j + 1] + a[j2 + 1]; - x1r = a[j] - a[j2]; - x1i = a[j + 1] - a[j2 + 1]; - y0r = a[j + 2] + a[j2 + 2]; - y0i = a[j + 3] + a[j2 + 3]; - y1r = a[j + 2] - a[j2 + 2]; - y1i = a[j + 3] - a[j2 + 3]; - x2r = a[j1] + a[j3]; - x2i = a[j1 + 1] + a[j3 + 1]; - x3r = a[j1] - a[j3]; - x3i = a[j1 + 1] - a[j3 + 1]; - y2r = a[j1 + 2] + a[j3 + 2]; - y2i = a[j1 + 3] + a[j3 + 3]; - y3r = a[j1 + 2] - a[j3 + 2]; - y3i = a[j1 + 3] - a[j3 + 3]; - a[j] = x0r + x2r; - a[j + 1] = x0i + x2i; - a[j + 2] = y0r + y2r; - a[j + 3] = y0i + y2i; - a[j1] = x0r - x2r; - a[j1 + 1] = x0i - x2i; - a[j1 + 2] = y0r - y2r; - a[j1 + 3] = y0i - y2i; - x0r = x1r - x3i; - x0i = x1i + x3r; - a[j2] = wk1r * x0r - wk1i * x0i; - a[j2 + 1] = wk1r * x0i + wk1i * x0r; - x0r = y1r - y3i; - x0i = y1i + y3r; - a[j2 + 2] = wd1r * x0r - wd1i * x0i; - a[j2 + 3] = wd1r * x0i + wd1i * x0r; - x0r = x1r + x3i; - x0i = x1i - x3r; - a[j3] = wk3r * x0r + wk3i * x0i; - a[j3 + 1] = wk3r * x0i - wk3i * x0r; - x0r = y1r + y3i; - x0i = y1i - y3r; - a[j3 + 2] = wd3r * x0r + wd3i * x0i; - a[j3 + 3] = wd3r * x0i - wd3i * x0r; - j0 = m - j; - j1 = j0 + m; - j2 = j1 + m; - j3 = j2 + m; - x0r = a[j0] + a[j2]; - x0i = a[j0 + 1] + a[j2 + 1]; - x1r = a[j0] - a[j2]; - x1i = a[j0 + 1] - a[j2 + 1]; - y0r = a[j0 - 2] + a[j2 - 2]; - y0i = a[j0 - 1] + a[j2 - 1]; - y1r = a[j0 - 2] - a[j2 - 2]; - y1i = a[j0 - 1] - a[j2 - 1]; - x2r = a[j1] + a[j3]; - x2i = a[j1 + 1] + a[j3 + 1]; - x3r = a[j1] - a[j3]; - x3i = a[j1 + 1] - a[j3 + 1]; - y2r = a[j1 - 2] + a[j3 - 2]; - y2i = a[j1 - 1] + a[j3 - 1]; - y3r = a[j1 - 2] - a[j3 - 2]; - y3i = a[j1 - 1] - a[j3 - 1]; - a[j0] = x0r + x2r; - a[j0 + 1] = x0i + x2i; - a[j0 - 2] = y0r + y2r; - a[j0 - 1] = y0i + y2i; - a[j1] = x0r - x2r; - a[j1 + 1] = x0i - x2i; - a[j1 - 2] = y0r - y2r; - a[j1 - 1] = y0i - y2i; - x0r = x1r - x3i; - x0i = x1i + x3r; - a[j2] = wk1i * x0r - wk1r * x0i; - a[j2 + 1] = wk1i * x0i + wk1r * x0r; - x0r = y1r - y3i; - x0i = y1i + y3r; - a[j2 - 2] = wd1i * x0r - wd1r * x0i; - a[j2 - 1] = wd1i * x0i + wd1r * x0r; - x0r = x1r + x3i; - x0i = x1i - x3r; - a[j3] = wk3i * x0r + wk3r * x0i; - a[j3 + 1] = wk3i * x0i - wk3r * x0r; - x0r = y1r + y3i; - x0i = y1i - y3r; - a[j3 - 2] = wd3i * x0r + wd3r * x0i; - a[j3 - 1] = wd3i * x0i - wd3r * x0r; - } - wk1r = csc1 * (wd1r + wn4r); - wk1i = csc1 * (wd1i + wn4r); - wk3r = csc3 * (wd3r - wn4r); - wk3i = csc3 * (wd3i - wn4r); - j0 = mh; - j1 = j0 + m; - j2 = j1 + m; - j3 = j2 + m; - x0r = a[j0 - 2] + a[j2 - 2]; - x0i = a[j0 - 1] + a[j2 - 1]; - x1r = a[j0 - 2] - a[j2 - 2]; - x1i = a[j0 - 1] - a[j2 - 1]; - x2r = a[j1 - 2] + a[j3 - 2]; - x2i = a[j1 - 1] + a[j3 - 1]; - x3r = a[j1 - 2] - a[j3 - 2]; - x3i = a[j1 - 1] - a[j3 - 1]; - a[j0 - 2] = x0r + x2r; - a[j0 - 1] = x0i + x2i; - a[j1 - 2] = x0r - x2r; - a[j1 - 1] = x0i - x2i; - x0r = x1r - x3i; - x0i = x1i + x3r; - a[j2 - 2] = wk1r * x0r - wk1i * x0i; - a[j2 - 1] = wk1r * x0i + wk1i * x0r; - x0r = x1r + x3i; - x0i = x1i - x3r; - a[j3 - 2] = wk3r * x0r + wk3i * x0i; - a[j3 - 1] = wk3r * x0i - wk3i * x0r; - x0r = a[j0] + a[j2]; - x0i = a[j0 + 1] + a[j2 + 1]; - x1r = a[j0] - a[j2]; - x1i = a[j0 + 1] - a[j2 + 1]; - x2r = a[j1] + a[j3]; - x2i = a[j1 + 1] + a[j3 + 1]; - x3r = a[j1] - a[j3]; - x3i = a[j1 + 1] - a[j3 + 1]; - a[j0] = x0r + x2r; - a[j0 + 1] = x0i + x2i; - a[j1] = x0r - x2r; - a[j1 + 1] = x0i - x2i; - x0r = x1r - x3i; - x0i = x1i + x3r; - a[j2] = wn4r * (x0r - x0i); - a[j2 + 1] = wn4r * (x0i + x0r); - x0r = x1r + x3i; - x0i = x1i - x3r; - a[j3] = -wn4r * (x0r + x0i); - a[j3 + 1] = -wn4r * (x0i - x0r); - x0r = a[j0 + 2] + a[j2 + 2]; - x0i = a[j0 + 3] + a[j2 + 3]; - x1r = a[j0 + 2] - a[j2 + 2]; - x1i = a[j0 + 3] - a[j2 + 3]; - x2r = a[j1 + 2] + a[j3 + 2]; - x2i = a[j1 + 3] + a[j3 + 3]; - x3r = a[j1 + 2] - a[j3 + 2]; - x3i = a[j1 + 3] - a[j3 + 3]; - a[j0 + 2] = x0r + x2r; - a[j0 + 3] = x0i + x2i; - a[j1 + 2] = x0r - x2r; - a[j1 + 3] = x0i - x2i; - x0r = x1r - x3i; - x0i = x1i + x3r; - a[j2 + 2] = wk1i * x0r - wk1r * x0i; - a[j2 + 3] = wk1i * x0i + wk1r * x0r; - x0r = x1r + x3i; - x0i = x1i - x3r; - a[j3 + 2] = wk3i * x0r + wk3r * x0i; - a[j3 + 3] = wk3i * x0i - wk3r * x0r; -} - - -void cftb1st(int n, double *a, double *w) -{ - int j, j0, j1, j2, j3, k, m, mh; - double wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i, - wd1r, wd1i, wd3r, wd3i; - double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, - y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i; - - mh = n >> 3; - m = 2 * mh; - j1 = m; - j2 = j1 + m; - j3 = j2 + m; - x0r = a[0] + a[j2]; - x0i = -a[1] - a[j2 + 1]; - x1r = a[0] - a[j2]; - x1i = -a[1] + a[j2 + 1]; - x2r = a[j1] + a[j3]; - x2i = a[j1 + 1] + a[j3 + 1]; - x3r = a[j1] - a[j3]; - x3i = a[j1 + 1] - a[j3 + 1]; - a[0] = x0r + x2r; - a[1] = x0i - x2i; - a[j1] = x0r - x2r; - a[j1 + 1] = x0i + x2i; - a[j2] = x1r + x3i; - a[j2 + 1] = x1i + x3r; - a[j3] = x1r - x3i; - a[j3 + 1] = x1i - x3r; - wn4r = w[1]; - csc1 = w[2]; - csc3 = w[3]; - wd1r = 1; - wd1i = 0; - wd3r = 1; - wd3i = 0; - k = 0; - for (j = 2; j < mh - 2; j += 4) { - k += 4; - wk1r = csc1 * (wd1r + w[k]); - wk1i = csc1 * (wd1i + w[k + 1]); - wk3r = csc3 * (wd3r + w[k + 2]); - wk3i = csc3 * (wd3i + w[k + 3]); - wd1r = w[k]; - wd1i = w[k + 1]; - wd3r = w[k + 2]; - wd3i = w[k + 3]; - j1 = j + m; - j2 = j1 + m; - j3 = j2 + m; - x0r = a[j] + a[j2]; - x0i = -a[j + 1] - a[j2 + 1]; - x1r = a[j] - a[j2]; - x1i = -a[j + 1] + a[j2 + 1]; - y0r = a[j + 2] + a[j2 + 2]; - y0i = -a[j + 3] - a[j2 + 3]; - y1r = a[j + 2] - a[j2 + 2]; - y1i = -a[j + 3] + a[j2 + 3]; - x2r = a[j1] + a[j3]; - x2i = a[j1 + 1] + a[j3 + 1]; - x3r = a[j1] - a[j3]; - x3i = a[j1 + 1] - a[j3 + 1]; - y2r = a[j1 + 2] + a[j3 + 2]; - y2i = a[j1 + 3] + a[j3 + 3]; - y3r = a[j1 + 2] - a[j3 + 2]; - y3i = a[j1 + 3] - a[j3 + 3]; - a[j] = x0r + x2r; - a[j + 1] = x0i - x2i; - a[j + 2] = y0r + y2r; - a[j + 3] = y0i - y2i; - a[j1] = x0r - x2r; - a[j1 + 1] = x0i + x2i; - a[j1 + 2] = y0r - y2r; - a[j1 + 3] = y0i + y2i; - x0r = x1r + x3i; - x0i = x1i + x3r; - a[j2] = wk1r * x0r - wk1i * x0i; - a[j2 + 1] = wk1r * x0i + wk1i * x0r; - x0r = y1r + y3i; - x0i = y1i + y3r; - a[j2 + 2] = wd1r * x0r - wd1i * x0i; - a[j2 + 3] = wd1r * x0i + wd1i * x0r; - x0r = x1r - x3i; - x0i = x1i - x3r; - a[j3] = wk3r * x0r + wk3i * x0i; - a[j3 + 1] = wk3r * x0i - wk3i * x0r; - x0r = y1r - y3i; - x0i = y1i - y3r; - a[j3 + 2] = wd3r * x0r + wd3i * x0i; - a[j3 + 3] = wd3r * x0i - wd3i * x0r; - j0 = m - j; - j1 = j0 + m; - j2 = j1 + m; - j3 = j2 + m; - x0r = a[j0] + a[j2]; - x0i = -a[j0 + 1] - a[j2 + 1]; - x1r = a[j0] - a[j2]; - x1i = -a[j0 + 1] + a[j2 + 1]; - y0r = a[j0 - 2] + a[j2 - 2]; - y0i = -a[j0 - 1] - a[j2 - 1]; - y1r = a[j0 - 2] - a[j2 - 2]; - y1i = -a[j0 - 1] + a[j2 - 1]; - x2r = a[j1] + a[j3]; - x2i = a[j1 + 1] + a[j3 + 1]; - x3r = a[j1] - a[j3]; - x3i = a[j1 + 1] - a[j3 + 1]; - y2r = a[j1 - 2] + a[j3 - 2]; - y2i = a[j1 - 1] + a[j3 - 1]; - y3r = a[j1 - 2] - a[j3 - 2]; - y3i = a[j1 - 1] - a[j3 - 1]; - a[j0] = x0r + x2r; - a[j0 + 1] = x0i - x2i; - a[j0 - 2] = y0r + y2r; - a[j0 - 1] = y0i - y2i; - a[j1] = x0r - x2r; - a[j1 + 1] = x0i + x2i; - a[j1 - 2] = y0r - y2r; - a[j1 - 1] = y0i + y2i; - x0r = x1r + x3i; - x0i = x1i + x3r; - a[j2] = wk1i * x0r - wk1r * x0i; - a[j2 + 1] = wk1i * x0i + wk1r * x0r; - x0r = y1r + y3i; - x0i = y1i + y3r; - a[j2 - 2] = wd1i * x0r - wd1r * x0i; - a[j2 - 1] = wd1i * x0i + wd1r * x0r; - x0r = x1r - x3i; - x0i = x1i - x3r; - a[j3] = wk3i * x0r + wk3r * x0i; - a[j3 + 1] = wk3i * x0i - wk3r * x0r; - x0r = y1r - y3i; - x0i = y1i - y3r; - a[j3 - 2] = wd3i * x0r + wd3r * x0i; - a[j3 - 1] = wd3i * x0i - wd3r * x0r; - } - wk1r = csc1 * (wd1r + wn4r); - wk1i = csc1 * (wd1i + wn4r); - wk3r = csc3 * (wd3r - wn4r); - wk3i = csc3 * (wd3i - wn4r); - j0 = mh; - j1 = j0 + m; - j2 = j1 + m; - j3 = j2 + m; - x0r = a[j0 - 2] + a[j2 - 2]; - x0i = -a[j0 - 1] - a[j2 - 1]; - x1r = a[j0 - 2] - a[j2 - 2]; - x1i = -a[j0 - 1] + a[j2 - 1]; - x2r = a[j1 - 2] + a[j3 - 2]; - x2i = a[j1 - 1] + a[j3 - 1]; - x3r = a[j1 - 2] - a[j3 - 2]; - x3i = a[j1 - 1] - a[j3 - 1]; - a[j0 - 2] = x0r + x2r; - a[j0 - 1] = x0i - x2i; - a[j1 - 2] = x0r - x2r; - a[j1 - 1] = x0i + x2i; - x0r = x1r + x3i; - x0i = x1i + x3r; - a[j2 - 2] = wk1r * x0r - wk1i * x0i; - a[j2 - 1] = wk1r * x0i + wk1i * x0r; - x0r = x1r - x3i; - x0i = x1i - x3r; - a[j3 - 2] = wk3r * x0r + wk3i * x0i; - a[j3 - 1] = wk3r * x0i - wk3i * x0r; - x0r = a[j0] + a[j2]; - x0i = -a[j0 + 1] - a[j2 + 1]; - x1r = a[j0] - a[j2]; - x1i = -a[j0 + 1] + a[j2 + 1]; - x2r = a[j1] + a[j3]; - x2i = a[j1 + 1] + a[j3 + 1]; - x3r = a[j1] - a[j3]; - x3i = a[j1 + 1] - a[j3 + 1]; - a[j0] = x0r + x2r; - a[j0 + 1] = x0i - x2i; - a[j1] = x0r - x2r; - a[j1 + 1] = x0i + x2i; - x0r = x1r + x3i; - x0i = x1i + x3r; - a[j2] = wn4r * (x0r - x0i); - a[j2 + 1] = wn4r * (x0i + x0r); - x0r = x1r - x3i; - x0i = x1i - x3r; - a[j3] = -wn4r * (x0r + x0i); - a[j3 + 1] = -wn4r * (x0i - x0r); - x0r = a[j0 + 2] + a[j2 + 2]; - x0i = -a[j0 + 3] - a[j2 + 3]; - x1r = a[j0 + 2] - a[j2 + 2]; - x1i = -a[j0 + 3] + a[j2 + 3]; - x2r = a[j1 + 2] + a[j3 + 2]; - x2i = a[j1 + 3] + a[j3 + 3]; - x3r = a[j1 + 2] - a[j3 + 2]; - x3i = a[j1 + 3] - a[j3 + 3]; - a[j0 + 2] = x0r + x2r; - a[j0 + 3] = x0i - x2i; - a[j1 + 2] = x0r - x2r; - a[j1 + 3] = x0i + x2i; - x0r = x1r + x3i; - x0i = x1i + x3r; - a[j2 + 2] = wk1i * x0r - wk1r * x0i; - a[j2 + 3] = wk1i * x0i + wk1r * x0r; - x0r = x1r - x3i; - x0i = x1i - x3r; - a[j3 + 2] = wk3i * x0r + wk3r * x0i; - a[j3 + 3] = wk3i * x0i - wk3r * x0r; -} - - -#ifdef USE_CDFT_THREADS -struct cdft_arg_st { - int n0; - int n; - double *a; - int nw; - double *w; -}; -typedef struct cdft_arg_st cdft_arg_t; - - -void cftrec4_th(int n, double *a, int nw, double *w) -{ - void *cftrec1_th(void *p); - void *cftrec2_th(void *p); - int i, idiv4, m, nthread; - cdft_thread_t th[4]; - cdft_arg_t ag[4]; - - nthread = 2; - idiv4 = 0; - m = n >> 1; - if (n > CDFT_4THREADS_BEGIN_N) { - nthread = 4; - idiv4 = 1; - m >>= 1; - } - for (i = 0; i < nthread; i++) { - ag[i].n0 = n; - ag[i].n = m; - ag[i].a = &a[i * m]; - ag[i].nw = nw; - ag[i].w = w; - if (i != idiv4) { - cdft_thread_create(&th[i], cftrec1_th, &ag[i]); - } else { - cdft_thread_create(&th[i], cftrec2_th, &ag[i]); - } - } - for (i = 0; i < nthread; i++) { - cdft_thread_wait(th[i]); - } -} - - -void *cftrec1_th(void *p) -{ - int cfttree(int n, int j, int k, double *a, int nw, double *w); - void cftleaf(int n, int isplt, double *a, int nw, double *w); - void cftmdl1(int n, double *a, double *w); - int isplt, j, k, m, n, n0, nw; - double *a, *w; - - n0 = ((cdft_arg_t *) p)->n0; - n = ((cdft_arg_t *) p)->n; - a = ((cdft_arg_t *) p)->a; - nw = ((cdft_arg_t *) p)->nw; - w = ((cdft_arg_t *) p)->w; - m = n0; - while (m > 512) { - m >>= 2; - cftmdl1(m, &a[n - m], &w[nw - (m >> 1)]); - } - cftleaf(m, 1, &a[n - m], nw, w); - k = 0; - for (j = n - m; j > 0; j -= m) { - k++; - isplt = cfttree(m, j, k, a, nw, w); - cftleaf(m, isplt, &a[j - m], nw, w); - } - return (void *) 0; -} - - -void *cftrec2_th(void *p) -{ - int cfttree(int n, int j, int k, double *a, int nw, double *w); - void cftleaf(int n, int isplt, double *a, int nw, double *w); - void cftmdl2(int n, double *a, double *w); - int isplt, j, k, m, n, n0, nw; - double *a, *w; - - n0 = ((cdft_arg_t *) p)->n0; - n = ((cdft_arg_t *) p)->n; - a = ((cdft_arg_t *) p)->a; - nw = ((cdft_arg_t *) p)->nw; - w = ((cdft_arg_t *) p)->w; - k = 1; - m = n0; - while (m > 512) { - m >>= 2; - k <<= 2; - cftmdl2(m, &a[n - m], &w[nw - m]); - } - cftleaf(m, 0, &a[n - m], nw, w); - k >>= 1; - for (j = n - m; j > 0; j -= m) { - k++; - isplt = cfttree(m, j, k, a, nw, w); - cftleaf(m, isplt, &a[j - m], nw, w); - } - return (void *) 0; -} -#endif /* USE_CDFT_THREADS */ - - -void cftrec4(int n, double *a, int nw, double *w) -{ - int cfttree(int n, int j, int k, double *a, int nw, double *w); - void cftleaf(int n, int isplt, double *a, int nw, double *w); - void cftmdl1(int n, double *a, double *w); - int isplt, j, k, m; - - m = n; - while (m > 512) { - m >>= 2; - cftmdl1(m, &a[n - m], &w[nw - (m >> 1)]); - } - cftleaf(m, 1, &a[n - m], nw, w); - k = 0; - for (j = n - m; j > 0; j -= m) { - k++; - isplt = cfttree(m, j, k, a, nw, w); - cftleaf(m, isplt, &a[j - m], nw, w); - } -} - - -int cfttree(int n, int j, int k, double *a, int nw, double *w) -{ - void cftmdl1(int n, double *a, double *w); - void cftmdl2(int n, double *a, double *w); - int i, isplt, m; - - if ((k & 3) != 0) { - isplt = k & 1; - if (isplt != 0) { - cftmdl1(n, &a[j - n], &w[nw - (n >> 1)]); - } else { - cftmdl2(n, &a[j - n], &w[nw - n]); - } - } else { - m = n; - for (i = k; (i & 3) == 0; i >>= 2) { - m <<= 2; - } - isplt = i & 1; - if (isplt != 0) { - while (m > 128) { - cftmdl1(m, &a[j - m], &w[nw - (m >> 1)]); - m >>= 2; - } - } else { - while (m > 128) { - cftmdl2(m, &a[j - m], &w[nw - m]); - m >>= 2; - } - } - } - return isplt; -} - - -void cftleaf(int n, int isplt, double *a, int nw, double *w) -{ - void cftmdl1(int n, double *a, double *w); - void cftmdl2(int n, double *a, double *w); - void cftf161(double *a, double *w); - void cftf162(double *a, double *w); - void cftf081(double *a, double *w); - void cftf082(double *a, double *w); - - if (n == 512) { - cftmdl1(128, a, &w[nw - 64]); - cftf161(a, &w[nw - 8]); - cftf162(&a[32], &w[nw - 32]); - cftf161(&a[64], &w[nw - 8]); - cftf161(&a[96], &w[nw - 8]); - cftmdl2(128, &a[128], &w[nw - 128]); - cftf161(&a[128], &w[nw - 8]); - cftf162(&a[160], &w[nw - 32]); - cftf161(&a[192], &w[nw - 8]); - cftf162(&a[224], &w[nw - 32]); - cftmdl1(128, &a[256], &w[nw - 64]); - cftf161(&a[256], &w[nw - 8]); - cftf162(&a[288], &w[nw - 32]); - cftf161(&a[320], &w[nw - 8]); - cftf161(&a[352], &w[nw - 8]); - if (isplt != 0) { - cftmdl1(128, &a[384], &w[nw - 64]); - cftf161(&a[480], &w[nw - 8]); - } else { - cftmdl2(128, &a[384], &w[nw - 128]); - cftf162(&a[480], &w[nw - 32]); - } - cftf161(&a[384], &w[nw - 8]); - cftf162(&a[416], &w[nw - 32]); - cftf161(&a[448], &w[nw - 8]); - } else { - cftmdl1(64, a, &w[nw - 32]); - cftf081(a, &w[nw - 8]); - cftf082(&a[16], &w[nw - 8]); - cftf081(&a[32], &w[nw - 8]); - cftf081(&a[48], &w[nw - 8]); - cftmdl2(64, &a[64], &w[nw - 64]); - cftf081(&a[64], &w[nw - 8]); - cftf082(&a[80], &w[nw - 8]); - cftf081(&a[96], &w[nw - 8]); - cftf082(&a[112], &w[nw - 8]); - cftmdl1(64, &a[128], &w[nw - 32]); - cftf081(&a[128], &w[nw - 8]); - cftf082(&a[144], &w[nw - 8]); - cftf081(&a[160], &w[nw - 8]); - cftf081(&a[176], &w[nw - 8]); - if (isplt != 0) { - cftmdl1(64, &a[192], &w[nw - 32]); - cftf081(&a[240], &w[nw - 8]); - } else { - cftmdl2(64, &a[192], &w[nw - 64]); - cftf082(&a[240], &w[nw - 8]); - } - cftf081(&a[192], &w[nw - 8]); - cftf082(&a[208], &w[nw - 8]); - cftf081(&a[224], &w[nw - 8]); - } -} - - -void cftmdl1(int n, double *a, double *w) -{ - int j, j0, j1, j2, j3, k, m, mh; - double wn4r, wk1r, wk1i, wk3r, wk3i; - double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; - - mh = n >> 3; - m = 2 * mh; - j1 = m; - j2 = j1 + m; - j3 = j2 + m; - x0r = a[0] + a[j2]; - x0i = a[1] + a[j2 + 1]; - x1r = a[0] - a[j2]; - x1i = a[1] - a[j2 + 1]; - x2r = a[j1] + a[j3]; - x2i = a[j1 + 1] + a[j3 + 1]; - x3r = a[j1] - a[j3]; - x3i = a[j1 + 1] - a[j3 + 1]; - a[0] = x0r + x2r; - a[1] = x0i + x2i; - a[j1] = x0r - x2r; - a[j1 + 1] = x0i - x2i; - a[j2] = x1r - x3i; - a[j2 + 1] = x1i + x3r; - a[j3] = x1r + x3i; - a[j3 + 1] = x1i - x3r; - wn4r = w[1]; - k = 0; - for (j = 2; j < mh; j += 2) { - k += 4; - wk1r = w[k]; - wk1i = w[k + 1]; - wk3r = w[k + 2]; - wk3i = w[k + 3]; - j1 = j + m; - j2 = j1 + m; - j3 = j2 + m; - x0r = a[j] + a[j2]; - x0i = a[j + 1] + a[j2 + 1]; - x1r = a[j] - a[j2]; - x1i = a[j + 1] - a[j2 + 1]; - x2r = a[j1] + a[j3]; - x2i = a[j1 + 1] + a[j3 + 1]; - x3r = a[j1] - a[j3]; - x3i = a[j1 + 1] - a[j3 + 1]; - a[j] = x0r + x2r; - a[j + 1] = x0i + x2i; - a[j1] = x0r - x2r; - a[j1 + 1] = x0i - x2i; - x0r = x1r - x3i; - x0i = x1i + x3r; - a[j2] = wk1r * x0r - wk1i * x0i; - a[j2 + 1] = wk1r * x0i + wk1i * x0r; - x0r = x1r + x3i; - x0i = x1i - x3r; - a[j3] = wk3r * x0r + wk3i * x0i; - a[j3 + 1] = wk3r * x0i - wk3i * x0r; - j0 = m - j; - j1 = j0 + m; - j2 = j1 + m; - j3 = j2 + m; - x0r = a[j0] + a[j2]; - x0i = a[j0 + 1] + a[j2 + 1]; - x1r = a[j0] - a[j2]; - x1i = a[j0 + 1] - a[j2 + 1]; - x2r = a[j1] + a[j3]; - x2i = a[j1 + 1] + a[j3 + 1]; - x3r = a[j1] - a[j3]; - x3i = a[j1 + 1] - a[j3 + 1]; - a[j0] = x0r + x2r; - a[j0 + 1] = x0i + x2i; - a[j1] = x0r - x2r; - a[j1 + 1] = x0i - x2i; - x0r = x1r - x3i; - x0i = x1i + x3r; - a[j2] = wk1i * x0r - wk1r * x0i; - a[j2 + 1] = wk1i * x0i + wk1r * x0r; - x0r = x1r + x3i; - x0i = x1i - x3r; - a[j3] = wk3i * x0r + wk3r * x0i; - a[j3 + 1] = wk3i * x0i - wk3r * x0r; - } - j0 = mh; - j1 = j0 + m; - j2 = j1 + m; - j3 = j2 + m; - x0r = a[j0] + a[j2]; - x0i = a[j0 + 1] + a[j2 + 1]; - x1r = a[j0] - a[j2]; - x1i = a[j0 + 1] - a[j2 + 1]; - x2r = a[j1] + a[j3]; - x2i = a[j1 + 1] + a[j3 + 1]; - x3r = a[j1] - a[j3]; - x3i = a[j1 + 1] - a[j3 + 1]; - a[j0] = x0r + x2r; - a[j0 + 1] = x0i + x2i; - a[j1] = x0r - x2r; - a[j1 + 1] = x0i - x2i; - x0r = x1r - x3i; - x0i = x1i + x3r; - a[j2] = wn4r * (x0r - x0i); - a[j2 + 1] = wn4r * (x0i + x0r); - x0r = x1r + x3i; - x0i = x1i - x3r; - a[j3] = -wn4r * (x0r + x0i); - a[j3 + 1] = -wn4r * (x0i - x0r); -} - - -void cftmdl2(int n, double *a, double *w) -{ - int j, j0, j1, j2, j3, k, kr, m, mh; - double wn4r, wk1r, wk1i, wk3r, wk3i, wd1r, wd1i, wd3r, wd3i; - double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, y0r, y0i, y2r, y2i; - - mh = n >> 3; - m = 2 * mh; - wn4r = w[1]; - j1 = m; - j2 = j1 + m; - j3 = j2 + m; - x0r = a[0] - a[j2 + 1]; - x0i = a[1] + a[j2]; - x1r = a[0] + a[j2 + 1]; - x1i = a[1] - a[j2]; - x2r = a[j1] - a[j3 + 1]; - x2i = a[j1 + 1] + a[j3]; - x3r = a[j1] + a[j3 + 1]; - x3i = a[j1 + 1] - a[j3]; - y0r = wn4r * (x2r - x2i); - y0i = wn4r * (x2i + x2r); - a[0] = x0r + y0r; - a[1] = x0i + y0i; - a[j1] = x0r - y0r; - a[j1 + 1] = x0i - y0i; - y0r = wn4r * (x3r - x3i); - y0i = wn4r * (x3i + x3r); - a[j2] = x1r - y0i; - a[j2 + 1] = x1i + y0r; - a[j3] = x1r + y0i; - a[j3 + 1] = x1i - y0r; - k = 0; - kr = 2 * m; - for (j = 2; j < mh; j += 2) { - k += 4; - wk1r = w[k]; - wk1i = w[k + 1]; - wk3r = w[k + 2]; - wk3i = w[k + 3]; - kr -= 4; - wd1i = w[kr]; - wd1r = w[kr + 1]; - wd3i = w[kr + 2]; - wd3r = w[kr + 3]; - j1 = j + m; - j2 = j1 + m; - j3 = j2 + m; - x0r = a[j] - a[j2 + 1]; - x0i = a[j + 1] + a[j2]; - x1r = a[j] + a[j2 + 1]; - x1i = a[j + 1] - a[j2]; - x2r = a[j1] - a[j3 + 1]; - x2i = a[j1 + 1] + a[j3]; - x3r = a[j1] + a[j3 + 1]; - x3i = a[j1 + 1] - a[j3]; - y0r = wk1r * x0r - wk1i * x0i; - y0i = wk1r * x0i + wk1i * x0r; - y2r = wd1r * x2r - wd1i * x2i; - y2i = wd1r * x2i + wd1i * x2r; - a[j] = y0r + y2r; - a[j + 1] = y0i + y2i; - a[j1] = y0r - y2r; - a[j1 + 1] = y0i - y2i; - y0r = wk3r * x1r + wk3i * x1i; - y0i = wk3r * x1i - wk3i * x1r; - y2r = wd3r * x3r + wd3i * x3i; - y2i = wd3r * x3i - wd3i * x3r; - a[j2] = y0r + y2r; - a[j2 + 1] = y0i + y2i; - a[j3] = y0r - y2r; - a[j3 + 1] = y0i - y2i; - j0 = m - j; - j1 = j0 + m; - j2 = j1 + m; - j3 = j2 + m; - x0r = a[j0] - a[j2 + 1]; - x0i = a[j0 + 1] + a[j2]; - x1r = a[j0] + a[j2 + 1]; - x1i = a[j0 + 1] - a[j2]; - x2r = a[j1] - a[j3 + 1]; - x2i = a[j1 + 1] + a[j3]; - x3r = a[j1] + a[j3 + 1]; - x3i = a[j1 + 1] - a[j3]; - y0r = wd1i * x0r - wd1r * x0i; - y0i = wd1i * x0i + wd1r * x0r; - y2r = wk1i * x2r - wk1r * x2i; - y2i = wk1i * x2i + wk1r * x2r; - a[j0] = y0r + y2r; - a[j0 + 1] = y0i + y2i; - a[j1] = y0r - y2r; - a[j1 + 1] = y0i - y2i; - y0r = wd3i * x1r + wd3r * x1i; - y0i = wd3i * x1i - wd3r * x1r; - y2r = wk3i * x3r + wk3r * x3i; - y2i = wk3i * x3i - wk3r * x3r; - a[j2] = y0r + y2r; - a[j2 + 1] = y0i + y2i; - a[j3] = y0r - y2r; - a[j3 + 1] = y0i - y2i; - } - wk1r = w[m]; - wk1i = w[m + 1]; - j0 = mh; - j1 = j0 + m; - j2 = j1 + m; - j3 = j2 + m; - x0r = a[j0] - a[j2 + 1]; - x0i = a[j0 + 1] + a[j2]; - x1r = a[j0] + a[j2 + 1]; - x1i = a[j0 + 1] - a[j2]; - x2r = a[j1] - a[j3 + 1]; - x2i = a[j1 + 1] + a[j3]; - x3r = a[j1] + a[j3 + 1]; - x3i = a[j1 + 1] - a[j3]; - y0r = wk1r * x0r - wk1i * x0i; - y0i = wk1r * x0i + wk1i * x0r; - y2r = wk1i * x2r - wk1r * x2i; - y2i = wk1i * x2i + wk1r * x2r; - a[j0] = y0r + y2r; - a[j0 + 1] = y0i + y2i; - a[j1] = y0r - y2r; - a[j1 + 1] = y0i - y2i; - y0r = wk1i * x1r - wk1r * x1i; - y0i = wk1i * x1i + wk1r * x1r; - y2r = wk1r * x3r - wk1i * x3i; - y2i = wk1r * x3i + wk1i * x3r; - a[j2] = y0r - y2r; - a[j2 + 1] = y0i - y2i; - a[j3] = y0r + y2r; - a[j3 + 1] = y0i + y2i; -} - - -void cftfx41(int n, double *a, int nw, double *w) -{ - void cftf161(double *a, double *w); - void cftf162(double *a, double *w); - void cftf081(double *a, double *w); - void cftf082(double *a, double *w); - - if (n == 128) { - cftf161(a, &w[nw - 8]); - cftf162(&a[32], &w[nw - 32]); - cftf161(&a[64], &w[nw - 8]); - cftf161(&a[96], &w[nw - 8]); - } else { - cftf081(a, &w[nw - 8]); - cftf082(&a[16], &w[nw - 8]); - cftf081(&a[32], &w[nw - 8]); - cftf081(&a[48], &w[nw - 8]); - } -} - - -void cftf161(double *a, double *w) -{ - double wn4r, wk1r, wk1i, - x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, - y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, - y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i, - y8r, y8i, y9r, y9i, y10r, y10i, y11r, y11i, - y12r, y12i, y13r, y13i, y14r, y14i, y15r, y15i; - - wn4r = w[1]; - wk1r = w[2]; - wk1i = w[3]; - x0r = a[0] + a[16]; - x0i = a[1] + a[17]; - x1r = a[0] - a[16]; - x1i = a[1] - a[17]; - x2r = a[8] + a[24]; - x2i = a[9] + a[25]; - x3r = a[8] - a[24]; - x3i = a[9] - a[25]; - y0r = x0r + x2r; - y0i = x0i + x2i; - y4r = x0r - x2r; - y4i = x0i - x2i; - y8r = x1r - x3i; - y8i = x1i + x3r; - y12r = x1r + x3i; - y12i = x1i - x3r; - x0r = a[2] + a[18]; - x0i = a[3] + a[19]; - x1r = a[2] - a[18]; - x1i = a[3] - a[19]; - x2r = a[10] + a[26]; - x2i = a[11] + a[27]; - x3r = a[10] - a[26]; - x3i = a[11] - a[27]; - y1r = x0r + x2r; - y1i = x0i + x2i; - y5r = x0r - x2r; - y5i = x0i - x2i; - x0r = x1r - x3i; - x0i = x1i + x3r; - y9r = wk1r * x0r - wk1i * x0i; - y9i = wk1r * x0i + wk1i * x0r; - x0r = x1r + x3i; - x0i = x1i - x3r; - y13r = wk1i * x0r - wk1r * x0i; - y13i = wk1i * x0i + wk1r * x0r; - x0r = a[4] + a[20]; - x0i = a[5] + a[21]; - x1r = a[4] - a[20]; - x1i = a[5] - a[21]; - x2r = a[12] + a[28]; - x2i = a[13] + a[29]; - x3r = a[12] - a[28]; - x3i = a[13] - a[29]; - y2r = x0r + x2r; - y2i = x0i + x2i; - y6r = x0r - x2r; - y6i = x0i - x2i; - x0r = x1r - x3i; - x0i = x1i + x3r; - y10r = wn4r * (x0r - x0i); - y10i = wn4r * (x0i + x0r); - x0r = x1r + x3i; - x0i = x1i - x3r; - y14r = wn4r * (x0r + x0i); - y14i = wn4r * (x0i - x0r); - x0r = a[6] + a[22]; - x0i = a[7] + a[23]; - x1r = a[6] - a[22]; - x1i = a[7] - a[23]; - x2r = a[14] + a[30]; - x2i = a[15] + a[31]; - x3r = a[14] - a[30]; - x3i = a[15] - a[31]; - y3r = x0r + x2r; - y3i = x0i + x2i; - y7r = x0r - x2r; - y7i = x0i - x2i; - x0r = x1r - x3i; - x0i = x1i + x3r; - y11r = wk1i * x0r - wk1r * x0i; - y11i = wk1i * x0i + wk1r * x0r; - x0r = x1r + x3i; - x0i = x1i - x3r; - y15r = wk1r * x0r - wk1i * x0i; - y15i = wk1r * x0i + wk1i * x0r; - x0r = y12r - y14r; - x0i = y12i - y14i; - x1r = y12r + y14r; - x1i = y12i + y14i; - x2r = y13r - y15r; - x2i = y13i - y15i; - x3r = y13r + y15r; - x3i = y13i + y15i; - a[24] = x0r + x2r; - a[25] = x0i + x2i; - a[26] = x0r - x2r; - a[27] = x0i - x2i; - a[28] = x1r - x3i; - a[29] = x1i + x3r; - a[30] = x1r + x3i; - a[31] = x1i - x3r; - x0r = y8r + y10r; - x0i = y8i + y10i; - x1r = y8r - y10r; - x1i = y8i - y10i; - x2r = y9r + y11r; - x2i = y9i + y11i; - x3r = y9r - y11r; - x3i = y9i - y11i; - a[16] = x0r + x2r; - a[17] = x0i + x2i; - a[18] = x0r - x2r; - a[19] = x0i - x2i; - a[20] = x1r - x3i; - a[21] = x1i + x3r; - a[22] = x1r + x3i; - a[23] = x1i - x3r; - x0r = y5r - y7i; - x0i = y5i + y7r; - x2r = wn4r * (x0r - x0i); - x2i = wn4r * (x0i + x0r); - x0r = y5r + y7i; - x0i = y5i - y7r; - x3r = wn4r * (x0r - x0i); - x3i = wn4r * (x0i + x0r); - x0r = y4r - y6i; - x0i = y4i + y6r; - x1r = y4r + y6i; - x1i = y4i - y6r; - a[8] = x0r + x2r; - a[9] = x0i + x2i; - a[10] = x0r - x2r; - a[11] = x0i - x2i; - a[12] = x1r - x3i; - a[13] = x1i + x3r; - a[14] = x1r + x3i; - a[15] = x1i - x3r; - x0r = y0r + y2r; - x0i = y0i + y2i; - x1r = y0r - y2r; - x1i = y0i - y2i; - x2r = y1r + y3r; - x2i = y1i + y3i; - x3r = y1r - y3r; - x3i = y1i - y3i; - a[0] = x0r + x2r; - a[1] = x0i + x2i; - a[2] = x0r - x2r; - a[3] = x0i - x2i; - a[4] = x1r - x3i; - a[5] = x1i + x3r; - a[6] = x1r + x3i; - a[7] = x1i - x3r; -} - - -void cftf162(double *a, double *w) -{ - double wn4r, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i, - x0r, x0i, x1r, x1i, x2r, x2i, - y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, - y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i, - y8r, y8i, y9r, y9i, y10r, y10i, y11r, y11i, - y12r, y12i, y13r, y13i, y14r, y14i, y15r, y15i; - - wn4r = w[1]; - wk1r = w[4]; - wk1i = w[5]; - wk3r = w[6]; - wk3i = -w[7]; - wk2r = w[8]; - wk2i = w[9]; - x1r = a[0] - a[17]; - x1i = a[1] + a[16]; - x0r = a[8] - a[25]; - x0i = a[9] + a[24]; - x2r = wn4r * (x0r - x0i); - x2i = wn4r * (x0i + x0r); - y0r = x1r + x2r; - y0i = x1i + x2i; - y4r = x1r - x2r; - y4i = x1i - x2i; - x1r = a[0] + a[17]; - x1i = a[1] - a[16]; - x0r = a[8] + a[25]; - x0i = a[9] - a[24]; - x2r = wn4r * (x0r - x0i); - x2i = wn4r * (x0i + x0r); - y8r = x1r - x2i; - y8i = x1i + x2r; - y12r = x1r + x2i; - y12i = x1i - x2r; - x0r = a[2] - a[19]; - x0i = a[3] + a[18]; - x1r = wk1r * x0r - wk1i * x0i; - x1i = wk1r * x0i + wk1i * x0r; - x0r = a[10] - a[27]; - x0i = a[11] + a[26]; - x2r = wk3i * x0r - wk3r * x0i; - x2i = wk3i * x0i + wk3r * x0r; - y1r = x1r + x2r; - y1i = x1i + x2i; - y5r = x1r - x2r; - y5i = x1i - x2i; - x0r = a[2] + a[19]; - x0i = a[3] - a[18]; - x1r = wk3r * x0r - wk3i * x0i; - x1i = wk3r * x0i + wk3i * x0r; - x0r = a[10] + a[27]; - x0i = a[11] - a[26]; - x2r = wk1r * x0r + wk1i * x0i; - x2i = wk1r * x0i - wk1i * x0r; - y9r = x1r - x2r; - y9i = x1i - x2i; - y13r = x1r + x2r; - y13i = x1i + x2i; - x0r = a[4] - a[21]; - x0i = a[5] + a[20]; - x1r = wk2r * x0r - wk2i * x0i; - x1i = wk2r * x0i + wk2i * x0r; - x0r = a[12] - a[29]; - x0i = a[13] + a[28]; - x2r = wk2i * x0r - wk2r * x0i; - x2i = wk2i * x0i + wk2r * x0r; - y2r = x1r + x2r; - y2i = x1i + x2i; - y6r = x1r - x2r; - y6i = x1i - x2i; - x0r = a[4] + a[21]; - x0i = a[5] - a[20]; - x1r = wk2i * x0r - wk2r * x0i; - x1i = wk2i * x0i + wk2r * x0r; - x0r = a[12] + a[29]; - x0i = a[13] - a[28]; - x2r = wk2r * x0r - wk2i * x0i; - x2i = wk2r * x0i + wk2i * x0r; - y10r = x1r - x2r; - y10i = x1i - x2i; - y14r = x1r + x2r; - y14i = x1i + x2i; - x0r = a[6] - a[23]; - x0i = a[7] + a[22]; - x1r = wk3r * x0r - wk3i * x0i; - x1i = wk3r * x0i + wk3i * x0r; - x0r = a[14] - a[31]; - x0i = a[15] + a[30]; - x2r = wk1i * x0r - wk1r * x0i; - x2i = wk1i * x0i + wk1r * x0r; - y3r = x1r + x2r; - y3i = x1i + x2i; - y7r = x1r - x2r; - y7i = x1i - x2i; - x0r = a[6] + a[23]; - x0i = a[7] - a[22]; - x1r = wk1i * x0r + wk1r * x0i; - x1i = wk1i * x0i - wk1r * x0r; - x0r = a[14] + a[31]; - x0i = a[15] - a[30]; - x2r = wk3i * x0r - wk3r * x0i; - x2i = wk3i * x0i + wk3r * x0r; - y11r = x1r + x2r; - y11i = x1i + x2i; - y15r = x1r - x2r; - y15i = x1i - x2i; - x1r = y0r + y2r; - x1i = y0i + y2i; - x2r = y1r + y3r; - x2i = y1i + y3i; - a[0] = x1r + x2r; - a[1] = x1i + x2i; - a[2] = x1r - x2r; - a[3] = x1i - x2i; - x1r = y0r - y2r; - x1i = y0i - y2i; - x2r = y1r - y3r; - x2i = y1i - y3i; - a[4] = x1r - x2i; - a[5] = x1i + x2r; - a[6] = x1r + x2i; - a[7] = x1i - x2r; - x1r = y4r - y6i; - x1i = y4i + y6r; - x0r = y5r - y7i; - x0i = y5i + y7r; - x2r = wn4r * (x0r - x0i); - x2i = wn4r * (x0i + x0r); - a[8] = x1r + x2r; - a[9] = x1i + x2i; - a[10] = x1r - x2r; - a[11] = x1i - x2i; - x1r = y4r + y6i; - x1i = y4i - y6r; - x0r = y5r + y7i; - x0i = y5i - y7r; - x2r = wn4r * (x0r - x0i); - x2i = wn4r * (x0i + x0r); - a[12] = x1r - x2i; - a[13] = x1i + x2r; - a[14] = x1r + x2i; - a[15] = x1i - x2r; - x1r = y8r + y10r; - x1i = y8i + y10i; - x2r = y9r - y11r; - x2i = y9i - y11i; - a[16] = x1r + x2r; - a[17] = x1i + x2i; - a[18] = x1r - x2r; - a[19] = x1i - x2i; - x1r = y8r - y10r; - x1i = y8i - y10i; - x2r = y9r + y11r; - x2i = y9i + y11i; - a[20] = x1r - x2i; - a[21] = x1i + x2r; - a[22] = x1r + x2i; - a[23] = x1i - x2r; - x1r = y12r - y14i; - x1i = y12i + y14r; - x0r = y13r + y15i; - x0i = y13i - y15r; - x2r = wn4r * (x0r - x0i); - x2i = wn4r * (x0i + x0r); - a[24] = x1r + x2r; - a[25] = x1i + x2i; - a[26] = x1r - x2r; - a[27] = x1i - x2i; - x1r = y12r + y14i; - x1i = y12i - y14r; - x0r = y13r - y15i; - x0i = y13i + y15r; - x2r = wn4r * (x0r - x0i); - x2i = wn4r * (x0i + x0r); - a[28] = x1r - x2i; - a[29] = x1i + x2r; - a[30] = x1r + x2i; - a[31] = x1i - x2r; -} - - -void cftf081(double *a, double *w) -{ - double wn4r, x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, - y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, - y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i; - - wn4r = w[1]; - x0r = a[0] + a[8]; - x0i = a[1] + a[9]; - x1r = a[0] - a[8]; - x1i = a[1] - a[9]; - x2r = a[4] + a[12]; - x2i = a[5] + a[13]; - x3r = a[4] - a[12]; - x3i = a[5] - a[13]; - y0r = x0r + x2r; - y0i = x0i + x2i; - y2r = x0r - x2r; - y2i = x0i - x2i; - y1r = x1r - x3i; - y1i = x1i + x3r; - y3r = x1r + x3i; - y3i = x1i - x3r; - x0r = a[2] + a[10]; - x0i = a[3] + a[11]; - x1r = a[2] - a[10]; - x1i = a[3] - a[11]; - x2r = a[6] + a[14]; - x2i = a[7] + a[15]; - x3r = a[6] - a[14]; - x3i = a[7] - a[15]; - y4r = x0r + x2r; - y4i = x0i + x2i; - y6r = x0r - x2r; - y6i = x0i - x2i; - x0r = x1r - x3i; - x0i = x1i + x3r; - x2r = x1r + x3i; - x2i = x1i - x3r; - y5r = wn4r * (x0r - x0i); - y5i = wn4r * (x0r + x0i); - y7r = wn4r * (x2r - x2i); - y7i = wn4r * (x2r + x2i); - a[8] = y1r + y5r; - a[9] = y1i + y5i; - a[10] = y1r - y5r; - a[11] = y1i - y5i; - a[12] = y3r - y7i; - a[13] = y3i + y7r; - a[14] = y3r + y7i; - a[15] = y3i - y7r; - a[0] = y0r + y4r; - a[1] = y0i + y4i; - a[2] = y0r - y4r; - a[3] = y0i - y4i; - a[4] = y2r - y6i; - a[5] = y2i + y6r; - a[6] = y2r + y6i; - a[7] = y2i - y6r; -} - - -void cftf082(double *a, double *w) -{ - double wn4r, wk1r, wk1i, x0r, x0i, x1r, x1i, - y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, - y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i; - - wn4r = w[1]; - wk1r = w[2]; - wk1i = w[3]; - y0r = a[0] - a[9]; - y0i = a[1] + a[8]; - y1r = a[0] + a[9]; - y1i = a[1] - a[8]; - x0r = a[4] - a[13]; - x0i = a[5] + a[12]; - y2r = wn4r * (x0r - x0i); - y2i = wn4r * (x0i + x0r); - x0r = a[4] + a[13]; - x0i = a[5] - a[12]; - y3r = wn4r * (x0r - x0i); - y3i = wn4r * (x0i + x0r); - x0r = a[2] - a[11]; - x0i = a[3] + a[10]; - y4r = wk1r * x0r - wk1i * x0i; - y4i = wk1r * x0i + wk1i * x0r; - x0r = a[2] + a[11]; - x0i = a[3] - a[10]; - y5r = wk1i * x0r - wk1r * x0i; - y5i = wk1i * x0i + wk1r * x0r; - x0r = a[6] - a[15]; - x0i = a[7] + a[14]; - y6r = wk1i * x0r - wk1r * x0i; - y6i = wk1i * x0i + wk1r * x0r; - x0r = a[6] + a[15]; - x0i = a[7] - a[14]; - y7r = wk1r * x0r - wk1i * x0i; - y7i = wk1r * x0i + wk1i * x0r; - x0r = y0r + y2r; - x0i = y0i + y2i; - x1r = y4r + y6r; - x1i = y4i + y6i; - a[0] = x0r + x1r; - a[1] = x0i + x1i; - a[2] = x0r - x1r; - a[3] = x0i - x1i; - x0r = y0r - y2r; - x0i = y0i - y2i; - x1r = y4r - y6r; - x1i = y4i - y6i; - a[4] = x0r - x1i; - a[5] = x0i + x1r; - a[6] = x0r + x1i; - a[7] = x0i - x1r; - x0r = y1r - y3i; - x0i = y1i + y3r; - x1r = y5r - y7r; - x1i = y5i - y7i; - a[8] = x0r + x1r; - a[9] = x0i + x1i; - a[10] = x0r - x1r; - a[11] = x0i - x1i; - x0r = y1r + y3i; - x0i = y1i - y3r; - x1r = y5r + y7r; - x1i = y5i + y7i; - a[12] = x0r - x1i; - a[13] = x0i + x1r; - a[14] = x0r + x1i; - a[15] = x0i - x1r; -} - - -void cftf040(double *a) -{ - double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; - - x0r = a[0] + a[4]; - x0i = a[1] + a[5]; - x1r = a[0] - a[4]; - x1i = a[1] - a[5]; - x2r = a[2] + a[6]; - x2i = a[3] + a[7]; - x3r = a[2] - a[6]; - x3i = a[3] - a[7]; - a[0] = x0r + x2r; - a[1] = x0i + x2i; - a[2] = x1r - x3i; - a[3] = x1i + x3r; - a[4] = x0r - x2r; - a[5] = x0i - x2i; - a[6] = x1r + x3i; - a[7] = x1i - x3r; -} - - -void cftb040(double *a) -{ - double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; - - x0r = a[0] + a[4]; - x0i = a[1] + a[5]; - x1r = a[0] - a[4]; - x1i = a[1] - a[5]; - x2r = a[2] + a[6]; - x2i = a[3] + a[7]; - x3r = a[2] - a[6]; - x3i = a[3] - a[7]; - a[0] = x0r + x2r; - a[1] = x0i + x2i; - a[2] = x1r + x3i; - a[3] = x1i - x3r; - a[4] = x0r - x2r; - a[5] = x0i - x2i; - a[6] = x1r - x3i; - a[7] = x1i + x3r; -} - - -void cftx020(double *a) -{ - double x0r, x0i; - - x0r = a[0] - a[2]; - x0i = a[1] - a[3]; - a[0] += a[2]; - a[1] += a[3]; - a[2] = x0r; - a[3] = x0i; -} - - -void rftfsub(int n, double *a, int nc, double *c) -{ - int j, k, kk, ks, m; - double wkr, wki, xr, xi, yr, yi; - - m = n >> 1; - ks = 2 * nc / m; - kk = 0; - for (j = 2; j < m; j += 2) { - k = n - j; - kk += ks; - wkr = 0.5 - c[nc - kk]; - wki = c[kk]; - xr = a[j] - a[k]; - xi = a[j + 1] + a[k + 1]; - yr = wkr * xr - wki * xi; - yi = wkr * xi + wki * xr; - a[j] -= yr; - a[j + 1] -= yi; - a[k] += yr; - a[k + 1] -= yi; - } -} - - -void rftbsub(int n, double *a, int nc, double *c) -{ - int j, k, kk, ks, m; - double wkr, wki, xr, xi, yr, yi; - - m = n >> 1; - ks = 2 * nc / m; - kk = 0; - for (j = 2; j < m; j += 2) { - k = n - j; - kk += ks; - wkr = 0.5 - c[nc - kk]; - wki = c[kk]; - xr = a[j] - a[k]; - xi = a[j + 1] + a[k + 1]; - yr = wkr * xr + wki * xi; - yi = wkr * xi - wki * xr; - a[j] -= yr; - a[j + 1] -= yi; - a[k] += yr; - a[k + 1] -= yi; - } -} - - -void dctsub(int n, double *a, int nc, double *c) -{ - int j, k, kk, ks, m; - double wkr, wki, xr; - - m = n >> 1; - ks = nc / n; - kk = 0; - for (j = 1; j < m; j++) { - k = n - j; - kk += ks; - wkr = c[kk] - c[nc - kk]; - wki = c[kk] + c[nc - kk]; - xr = wki * a[j] - wkr * a[k]; - a[j] = wkr * a[j] + wki * a[k]; - a[k] = xr; - } - a[m] *= c[0]; -} - - -void dstsub(int n, double *a, int nc, double *c) -{ - int j, k, kk, ks, m; - double wkr, wki, xr; - - m = n >> 1; - ks = nc / n; - kk = 0; - for (j = 1; j < m; j++) { - k = n - j; - kk += ks; - wkr = c[kk] - c[nc - kk]; - wki = c[kk] + c[nc - kk]; - xr = wki * a[k] - wkr * a[j]; - a[k] = wkr * a[k] + wki * a[j]; - a[j] = xr; - } - a[m] *= c[0]; -} - diff --git a/src/libprojectM/Audio/fftsg.h b/src/libprojectM/Audio/fftsg.h deleted file mode 100755 index a83669cc23..0000000000 --- a/src/libprojectM/Audio/fftsg.h +++ /dev/null @@ -1,35 +0,0 @@ -/** - * projectM -- Milkdrop-esque visualisation SDK - * Copyright (C)2003-2007 projectM Team - * - * This library is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * This library is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with this library; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - * See 'LICENSE.txt' included within this release - * - */ -/** - * $Id: fftsg.h,v 1.1.1.1 2005/12/23 18:05:00 psperl Exp $ - * - * Wrapper for rdft() and friends - * - * $Log$ - */ - -#ifndef _FFTSG_H -#define _FFTSG_H - -extern void rdft(int n, int isgn, double *a, int *ip, double *w); - -#endif /** !_FFTSG_H */ - diff --git a/src/libprojectM/CMakeLists.txt b/src/libprojectM/CMakeLists.txt index 8225212c49..ee502e5160 100644 --- a/src/libprojectM/CMakeLists.txt +++ b/src/libprojectM/CMakeLists.txt @@ -12,20 +12,12 @@ if(CMAKE_SYSTEM_NAME STREQUAL "Windows") include_directories("${MSVC_EXTRA_INCLUDE_DIR}") endif() +add_subdirectory(Audio) add_subdirectory(MilkdropPreset) add_subdirectory(Renderer) add_library(projectM_main OBJECT "${PROJECTM_EXPORT_HEADER}" - Audio/AudioConstants.hpp - Audio/BeatDetect.cpp - Audio/BeatDetect.hpp - Audio/FrameAudioData.cpp - Audio/FrameAudioData.hpp - Audio/PCM.cpp - Audio/PCM.hpp - Audio/fftsg.cpp - Audio/fftsg.h Preset.hpp PresetFactory.cpp PresetFactory.hpp @@ -43,6 +35,7 @@ add_library(projectM_main OBJECT target_link_libraries(projectM_main PUBLIC + Audio MilkdropPreset Renderer hlslparser @@ -71,6 +64,7 @@ target_include_directories(projectM_main # This syntax will pull in the compiled object files into the final library. add_library(projectM ${PROJECTM_DUMMY_SOURCE_FILE} # CMake needs at least one "real" source file. + $ $ $ $ diff --git a/src/libprojectM/MilkdropPreset/Constants.hpp b/src/libprojectM/MilkdropPreset/Constants.hpp index 0229601f1e..1d9d71cd94 100644 --- a/src/libprojectM/MilkdropPreset/Constants.hpp +++ b/src/libprojectM/MilkdropPreset/Constants.hpp @@ -12,5 +12,4 @@ static constexpr int TVarCount = 8; //!< Number of T variables available. static constexpr int CustomWaveformCount = 4; //!< Number of custom waveforms (expression-driven) which can be used in a preset. static constexpr int CustomShapeCount = 4; //!< Number of custom shapes (expression-driven) which can be used in a preset. -static constexpr int RenderWaveformSamples = 480; //!< Number of custom waveform data samples available for rendering a frame. static constexpr int WaveformMaxPoints = 512; //!< Maximum number of waveform points. diff --git a/src/libprojectM/MilkdropPreset/CustomShape.cpp b/src/libprojectM/MilkdropPreset/CustomShape.cpp index 5d40b26644..7a3405f898 100644 --- a/src/libprojectM/MilkdropPreset/CustomShape.cpp +++ b/src/libprojectM/MilkdropPreset/CustomShape.cpp @@ -2,7 +2,6 @@ #include "PresetFileParser.hpp" -#include "Audio/BeatDetect.hpp" #include "Renderer/TextureManager.hpp" #include diff --git a/src/libprojectM/MilkdropPreset/CustomWaveform.cpp b/src/libprojectM/MilkdropPreset/CustomWaveform.cpp index c34d09638d..2463bfddd6 100644 --- a/src/libprojectM/MilkdropPreset/CustomWaveform.cpp +++ b/src/libprojectM/MilkdropPreset/CustomWaveform.cpp @@ -12,7 +12,7 @@ #include #include -static constexpr int CustomWaveformMaxSamples = std::max(RenderWaveformSamples, libprojectM::Audio::SpectrumSamples); +static constexpr int CustomWaveformMaxSamples = std::max(libprojectM::Audio::WaveformSamples, libprojectM::Audio::SpectrumSamples); CustomWaveform::CustomWaveform(PresetState& presetState) : RenderItem() @@ -73,17 +73,15 @@ void CustomWaveform::CompileCodeAndRunInitExpressions(const PerFrameContext& pre void CustomWaveform::Draw(const PerFrameContext& presetPerFrameContext) { - // Some safety assertions if someone plays around and changes the values in the wrong way. - static_assert(WaveformMaxPoints <= libprojectM::Audio::SpectrumSamples, "CustomWaveformMaxPoints is larger than SpectrumSamples"); - static_assert(WaveformMaxPoints <= libprojectM::Audio::WaveformSamples, "CustomWaveformMaxPoints is larger than WaveformSamples"); - static_assert(RenderWaveformSamples <= WaveformMaxPoints, "CustomWaveformSamples is larger than CustomWaveformMaxPoints"); + static_assert(libprojectM::Audio::WaveformSamples <= WaveformMaxPoints, "WaveformMaxPoints is larger than WaveformSamples"); + static_assert(libprojectM::Audio::SpectrumSamples <= WaveformMaxPoints, "WaveformMaxPoints is larger than SpectrumSamples"); if (!m_enabled) { return; } - int const maxSampleCount{m_spectrum ? libprojectM::Audio::SpectrumSamples : RenderWaveformSamples}; + int const maxSampleCount{m_spectrum ? libprojectM::Audio::SpectrumSamples : libprojectM::Audio::WaveformSamples}; // Initialize and execute per-frame code LoadPerFrameEvaluationVariables(presetPerFrameContext); @@ -92,7 +90,8 @@ void CustomWaveform::Draw(const PerFrameContext& presetPerFrameContext) // Copy Q and T vars to per-point context InitPerPointEvaluationVariables(); - int sampleCount = std::min(WaveformMaxPoints, static_cast(*m_perFrameContext.samples)); + int sampleCount = std::min(maxSampleCount, static_cast(*m_perFrameContext.samples)); + sampleCount -= m_sep; // If there aren't enough samples to draw a single line or dot, skip drawing the waveform. if ((m_useDots && sampleCount < 1) || sampleCount < 2) @@ -107,8 +106,8 @@ void CustomWaveform::Draw(const PerFrameContext& presetPerFrameContext) ? m_presetState.audioData.spectrumRight.data() : m_presetState.audioData.waveformRight.data(); - //const float mult = m_scaling * m_presetState.waveScale * (m_spectrum ? 0.15f : 0.004f); - const float mult = m_scaling * m_presetState.waveScale * (m_spectrum ? 0.05f : 1.0f); + const float mult = m_scaling * m_presetState.waveScale * (m_spectrum ? 0.15f : 0.004f); + //const float mult = m_scaling * m_presetState.waveScale * (m_spectrum ? 0.05f : 1.0f); // PCM data smoothing const int offset1 = m_spectrum ? 0 : (maxSampleCount - sampleCount) / 2 - m_sep / 2; diff --git a/src/libprojectM/MilkdropPreset/Waveform.cpp b/src/libprojectM/MilkdropPreset/Waveform.cpp index 963672f5b2..eedea82663 100644 --- a/src/libprojectM/MilkdropPreset/Waveform.cpp +++ b/src/libprojectM/MilkdropPreset/Waveform.cpp @@ -3,9 +3,10 @@ #include "PerFrameContext.hpp" #include "PresetState.hpp" -#include "Audio/BeatDetect.hpp" #include "projectM-opengl.h" +#include "Audio/AudioConstants.hpp" + #include #include @@ -252,7 +253,11 @@ void Waveform::MaximizeColors(const PerFrameContext& presetPerFrameContext) void Waveform::WaveformMath(const PerFrameContext& presetPerFrameContext) { - constexpr size_t audioSamples{512}; + static_assert(libprojectM::Audio::WaveformSamples <= WaveformMaxPoints, "WaveformMaxPoints is larger than WaveformSamples"); + + using libprojectM::Audio::WaveformSamples; + + constexpr size_t audioSamples{WaveformMaxPoints}; delete[] m_wave1Vertices; m_wave1Vertices = nullptr; @@ -284,15 +289,16 @@ void Waveform::WaveformMath(const PerFrameContext& presetPerFrameContext) // Tie size of waveform to beatSensitivity // ToDo: Beat sensitivity was probably not the correct value here? - const float volumeScale = m_presetState.audioData.vol; - if (volumeScale != 1.0) + //const float volumeScale = m_presetState.audioData.vol; + //if (volumeScale != 1.0) + //{ + // Scale waveform data to -1 .. 1 range + for (size_t i = 0; i < pcmDataL.size(); ++i) { - for (size_t i = 0; i < pcmDataL.size(); ++i) - { - pcmDataL[i] *= volumeScale; - pcmDataR[i] *= volumeScale; - } + pcmDataL[i] /= 128.0f; + pcmDataR[i] /= 128.0f; } + //} // Aspect multipliers float aspectX{1.0f}; @@ -328,11 +334,11 @@ void Waveform::WaveformMath(const PerFrameContext& presetPerFrameContext) case Mode::Circle: { m_loop = true; - m_samples = RenderWaveformSamples / 2; + m_samples = WaveformSamples / 2; m_wave1Vertices = new Point[m_samples](); - const int sampleOffset{(RenderWaveformSamples - m_samples) / 2}; + const int sampleOffset{(WaveformSamples - m_samples) / 2}; const float inverseSamplesMinusOne{1.0f / static_cast(m_samples)}; @@ -356,7 +362,7 @@ void Waveform::WaveformMath(const PerFrameContext& presetPerFrameContext) case Mode::XYOscillationSpiral: //circularly moving waveform { - m_samples = RenderWaveformSamples / 2; + m_samples = WaveformSamples / 2; m_wave1Vertices = new Point[m_samples](); @@ -374,7 +380,7 @@ void Waveform::WaveformMath(const PerFrameContext& presetPerFrameContext) case Mode::CenteredSpiro: case Mode::CenteredSpiroVolume: { // Both "centered spiro" waveforms are identical. Only difference is the alpha value. // Alpha calculation is handled in MaximizeColors(). - m_samples = RenderWaveformSamples; + m_samples = WaveformSamples; m_wave1Vertices = new Point[m_samples](); @@ -388,7 +394,7 @@ void Waveform::WaveformMath(const PerFrameContext& presetPerFrameContext) } case Mode::DerivativeLine: { - m_samples = RenderWaveformSamples; + m_samples = WaveformSamples; if (m_samples > m_presetState.renderContext.viewportSizeX / 3) { @@ -397,7 +403,7 @@ void Waveform::WaveformMath(const PerFrameContext& presetPerFrameContext) m_wave1Vertices = new Point[m_samples](); - int const sampleOffset = (RenderWaveformSamples - m_samples) / 2; + int const sampleOffset = (WaveformSamples - m_samples) / 2; const float w1 = 0.45f + 0.5f * (mysteryWaveParam * 0.5f + 0.5f); const float w2 = 1.0f - w1; @@ -422,7 +428,7 @@ void Waveform::WaveformMath(const PerFrameContext& presetPerFrameContext) } case Mode::ExplosiveHash: { - m_samples = RenderWaveformSamples; + m_samples = WaveformSamples; m_wave1Vertices = new Point[m_samples](); @@ -449,7 +455,7 @@ void Waveform::WaveformMath(const PerFrameContext& presetPerFrameContext) } else { - m_samples = RenderWaveformSamples / 2; + m_samples = WaveformSamples / 2; if (m_samples > m_presetState.renderContext.viewportSizeX / 3) { @@ -463,7 +469,7 @@ void Waveform::WaveformMath(const PerFrameContext& presetPerFrameContext) m_wave2Vertices = new Point[m_samples](); } - const int sampleOffset = (RenderWaveformSamples - m_samples) / 2; + const int sampleOffset = (WaveformSamples - m_samples) / 2; const float angle = 1.57f * mysteryWaveParam; // from -PI/2 to PI/2 float dx = cosf(angle); diff --git a/src/libprojectM/ProjectM.cpp b/src/libprojectM/ProjectM.cpp index 9400cc73f5..799a99639f 100644 --- a/src/libprojectM/ProjectM.cpp +++ b/src/libprojectM/ProjectM.cpp @@ -25,7 +25,6 @@ #include "PresetFactoryManager.hpp" #include "TimeKeeper.hpp" -#include "Audio/BeatDetect.hpp" #include "Audio/PCM.hpp" //Sound data handler (buffering, FFT, etc.) #include "Renderer/CopyTexture.hpp" @@ -143,13 +142,16 @@ void ProjectM::RenderFrame() return; } - m_timeKeeper->UpdateTimers(); - m_beatDetect->CalculateBeatStatistics(); - #if PROJECTM_USE_THREADS std::lock_guard guard(m_presetSwitchMutex); #endif + // Update FPS and other timer values. + m_timeKeeper->UpdateTimers(); + + // Update and retrieve audio data + m_audioStorage.UpdateFrameAudioData(m_timeKeeper->SecondsSinceLastFrame(), m_frameCount); + auto audioData = m_audioStorage.GetFrameAudioData(); // Check if the preset isn't locked, and we've not already notified the user if (!m_presetChangeNotified) @@ -161,7 +163,8 @@ void ProjectM::RenderFrame() PresetSwitchRequestedEvent(false); } else if (m_hardCutEnabled && - (m_beatDetect->vol - m_beatDetect->volOld > m_hardCutSensitivity) && + m_frameCount > 50 && + (audioData.vol - m_previousFrameVolume > m_hardCutSensitivity) && m_timeKeeper->CanHardCut()) { m_presetChangeNotified = true; @@ -198,7 +201,6 @@ void ProjectM::RenderFrame() } auto renderContext = GetRenderContext(); - auto audioData = m_beatDetect->GetFrameAudioData(); if (m_transition != nullptr && m_transitioningPreset != nullptr) { @@ -231,6 +233,7 @@ void ProjectM::RenderFrame() } m_frameCount++; + m_previousFrameVolume = audioData.vol; } void ProjectM::Initialize() @@ -245,10 +248,6 @@ void ProjectM::Initialize() /** Initialise per-pixel matrix calculations */ /** We need to initialise this before the builtin param db otherwise bass/mid etc won't bind correctly */ - assert(!m_beatDetect); - - m_beatDetect = std::make_unique(m_pcm); - m_textureManager = std::make_unique(m_textureSearchPaths); m_transitionShaderManager = std::make_unique(); @@ -260,7 +259,6 @@ void ProjectM::Initialize() /* Set the seed to the current time in seconds */ srand(time(nullptr)); - ResetEngine(); LoadIdlePreset(); #if PROJECTM_USE_THREADS @@ -277,17 +275,6 @@ void ProjectM::LoadIdlePreset() assert(m_activePreset); } -/* Reinitializes the engine variables to a default (conservative and sane) value */ -void ProjectM::ResetEngine() -{ - - if (m_beatDetect != NULL) - { - m_beatDetect->Reset(); - m_beatDetect->beatSensitivity = m_beatSensitivity; - } -} - /** Resets OpenGL state */ void ProjectM::ResetOpenGL(size_t width, size_t height) { @@ -361,12 +348,12 @@ auto ProjectM::PresetLocked() const -> bool void ProjectM::SetBeatSensitivity(float sensitivity) { - m_beatDetect->beatSensitivity = std::min(std::max(0.0f, sensitivity), 2.0f); + m_beatSensitivity = std::min(std::max(0.0f, sensitivity), 2.0f); } auto ProjectM::GetBeatSensitivity() const -> float { - return m_beatDetect->beatSensitivity; + return m_beatSensitivity; } auto ProjectM::SoftCutDuration() const -> double @@ -481,7 +468,7 @@ void ProjectM::SetMeshSize(size_t meshResolutionX, size_t meshResolutionY) auto ProjectM::PCM() -> libprojectM::Audio::PCM& { - return m_pcm; + return m_audioStorage; } void ProjectM::Touch(float touchX, float touchY, int pressure, int touchType) diff --git a/src/libprojectM/ProjectM.hpp b/src/libprojectM/ProjectM.hpp index f850b33643..bf5cae969b 100644 --- a/src/libprojectM/ProjectM.hpp +++ b/src/libprojectM/ProjectM.hpp @@ -212,8 +212,6 @@ class PROJECTM_EXPORT ProjectM private: void Initialize(); - void ResetEngine(); - void StartPresetTransition(std::unique_ptr&& preset, bool hardCut); void LoadIdlePreset(); @@ -226,7 +224,6 @@ class PROJECTM_EXPORT ProjectM #endif - class libprojectM::Audio::PCM m_pcm; //!< Audio data buffer and analyzer instance. size_t m_meshX{32}; //!< Per-point mesh horizontal resolution. size_t m_meshY{24}; //!< Per-point mesh vertical resolution. @@ -241,6 +238,7 @@ class PROJECTM_EXPORT ProjectM float m_beatSensitivity{1.0}; //!< General beat sensitivity modifier for presets. bool m_aspectCorrection{true}; //!< If true, corrects aspect ratio for non-rectangular windows. float m_easterEgg{0.0}; //!< Random preset duration modifier. See TimeKeeper class. + float m_previousFrameVolume{}; //!< Volume in previous frame, used for hard cuts. std::vector m_textureSearchPaths; ///!< List of paths to search for texture files @@ -252,10 +250,10 @@ class PROJECTM_EXPORT ProjectM std::unique_ptr m_presetFactoryManager; //!< Provides access to all available preset factories. + libprojectM::Audio::PCM m_audioStorage; //!< Audio data buffer and analyzer instance. std::unique_ptr m_textureManager; //!< The texture manager. std::unique_ptr m_transitionShaderManager; //!< The transition shader manager. std::unique_ptr m_textureCopier; //!< Class that copies textures 1:1 to another texture or framebuffer. - std::unique_ptr m_beatDetect; //!< The beat detection class. std::unique_ptr m_activePreset; //!< Currently loaded preset. std::unique_ptr m_transitioningPreset; //!< Destination preset when smooth preset switching. std::unique_ptr m_transition; //!< Transition effect used for blending. diff --git a/src/libprojectM/ProjectMCWrapper.cpp b/src/libprojectM/ProjectMCWrapper.cpp index 7e5da6c299..3baabe5fd6 100644 --- a/src/libprojectM/ProjectMCWrapper.cpp +++ b/src/libprojectM/ProjectMCWrapper.cpp @@ -2,8 +2,9 @@ #include "projectM-4/projectM.h" +#include "Audio/AudioConstants.hpp" + #include -#include #include void projectMWrapper::PresetSwitchRequestedEvent(bool isHardCut) const @@ -337,7 +338,7 @@ void projectm_set_window_size(projectm_handle instance, size_t width, size_t hei unsigned int projectm_pcm_get_max_samples() { - return libprojectM::Audio::PCM::maxSamples; + return libprojectM::Audio::WaveformSamples; } template @@ -345,14 +346,7 @@ static auto PcmAdd(projectm_handle instance, const BufferType* samples, unsigned { auto* projectMInstance = handle_to_instance(instance); - if (channels == PROJECTM_MONO) - { - projectMInstance->PCM().AddMono(samples, count); - } - else - { - projectMInstance->PCM().AddStereo(samples, count); - } + projectMInstance->PCM().Add(samples, channels, count); } auto projectm_pcm_add_float(projectm_handle instance, const float* samples, unsigned int count, projectm_channels channels) -> void diff --git a/src/libprojectM/TimeKeeper.cpp b/src/libprojectM/TimeKeeper.cpp index 839401e6c7..91c35a0edf 100644 --- a/src/libprojectM/TimeKeeper.cpp +++ b/src/libprojectM/TimeKeeper.cpp @@ -17,7 +17,9 @@ void TimeKeeper::UpdateTimers() { auto currentTime = std::chrono::high_resolution_clock::now(); - m_currentTime = std::chrono::duration(currentTime - m_startTime).count(); + double currentFrameTime = std::chrono::duration(currentTime - m_startTime).count(); + m_secondsSinceLastFrame = currentFrameTime - m_currentTime; + m_currentTime = currentFrameTime; m_presetFrameA++; m_presetFrameB++; } diff --git a/src/libprojectM/TimeKeeper.hpp b/src/libprojectM/TimeKeeper.hpp index c79e6c4884..273f1c6ec7 100644 --- a/src/libprojectM/TimeKeeper.hpp +++ b/src/libprojectM/TimeKeeper.hpp @@ -6,7 +6,6 @@ class TimeKeeper { public: - TimeKeeper(double presetDuration, double smoothDuration, double hardcutDuration, double easterEgg); void UpdateTimers(); @@ -79,27 +78,31 @@ class TimeKeeper m_easterEgg = value; } -private: + inline auto SecondsSinceLastFrame() const -> double + { + return m_secondsSinceLastFrame; + } +private: /* The first ticks value of the application */ - std::chrono::high_resolution_clock::time_point m_startTime{ std::chrono::high_resolution_clock::now() }; - - double m_easterEgg{ 0.0 }; + std::chrono::high_resolution_clock::time_point m_startTime{std::chrono::high_resolution_clock::now()}; - double m_presetDuration{ 0.0 }; - double m_presetDurationA{ 0.0 }; - double m_presetDurationB{ 0.0 }; - double m_softCutDuration{ 0.0 }; - double m_hardCutDuration{ 0.0 }; + double m_secondsSinceLastFrame{}; - double m_currentTime{ 0.0 }; - double m_presetTimeA{ 0.0 }; - double m_presetTimeB{ 0.0 }; + double m_easterEgg{}; - int m_presetFrameA{ 0 }; - int m_presetFrameB{ 0 }; + double m_presetDuration{}; + double m_presetDurationA{}; + double m_presetDurationB{}; + double m_softCutDuration{}; + double m_hardCutDuration{}; - bool m_isSmoothing{ false }; + double m_currentTime{}; + double m_presetTimeA{}; + double m_presetTimeB{}; + int m_presetFrameA{}; + int m_presetFrameB{}; + bool m_isSmoothing{false}; }; diff --git a/tests/libprojectM/CMakeLists.txt b/tests/libprojectM/CMakeLists.txt index 4dc90eae82..963818ee24 100644 --- a/tests/libprojectM/CMakeLists.txt +++ b/tests/libprojectM/CMakeLists.txt @@ -1,9 +1,9 @@ find_package(GTest 1.10 REQUIRED NO_MODULE) add_executable(projectM-unittest - PCMTest.cpp PresetFileParserTest.cpp + $ $ $ $ diff --git a/tests/libprojectM/PCMTest.cpp b/tests/libprojectM/PCMTest.cpp deleted file mode 100644 index 7a7e83dd94..0000000000 --- a/tests/libprojectM/PCMTest.cpp +++ /dev/null @@ -1,140 +0,0 @@ -#include - -#include