見出し画像

Copilot君と一緒に頑張ったリサンプリングプログラム(完)

この世界にYouTubeが流行している限り、44.1KHz・48KHzの両対応から逃れられない。
ちょっといいDAC買ったのにお風呂サウンドになっちゃうのも、音が刺さるのも、高音寄りな気がするのも、リサンプリングのせいだったりする。
リサンプリングの限界を知って、システムの方向を見直した方が幸せになれるかも?
世の中のリサンプラーはSinc Lenが512くらいで高品質らしいけど、計算したらだいたいSinc Len 80万くらい必要そうだったので、全力のやつ作ってみた。

2025/1/18 22:05
計算ミスがあったので修正しました。
2025/1/19 5:19
ローパスモード(-lp)を追加しました。
2025/1/19 8:58
ローパスモード(-lp)の実装を変更しました。
リサンプリングの実装を変更し、元のサンプルの利用せず、全サンプルをSinc関数を通して計算するようにしました。


機能

  • 音声ファイルをリサンプリングします

  • Sinc関数で全サンプルの再計算

  • Sinc関数を利用したローパスフィルター

  • 窓パラメーターの自動計算

    1. カイザー窓のβをサイドローブが出力ビット深度以下になるように計算

    2. カイザー窓の誤差精度を1Hz未満に収まるようにSinc Lenを計算

    3. カイザー窓とSinc関数の窓長を2 * Sinc Len + 1に設定

  • 常に1.0秒単位になるように無音を追加してファイル長を統一

  • 窓関数とSinc関数の乗算時点では80-bit浮動小数点の精度
    その後の音声信号との乗算は64-bit浮動小数点の精度

  • CUDAを使って爆速
    整数倍 x8オーバーサンプリングで約1時間(RTX3060)
    非整数倍 x8.7オーバーサンプリングで約1日(RTX3060)


動作環境

i5-13500 / RTX3060 / Debian

Windowsの場合は無理せずWSLにUbuntuを入れれば同じようにできるはず。
その際はWSL専用のNvidiaドライバーがあるっぽい。
途中でエラーになったらCopilotに聞けばだいたいなんとかなる。
CPUはAMDでも動くんじゃないかな。


インストールしとくもの

nvidia関係のやつでひっかかったらCopilot君に聞いてください

sudo apt update && \
sudo apt install libsndfile1 fftw3 && \
sudo apt install nvidia-driver firmware-misc-nonfree nvidia-cuda-toolkit nvidia-cuda-dev nvidia-cudnn

ファイルの準備

記事の最後にコードを載せてるのでそれを手作業で保存。
中身見えなくても怖くない人はダウンロード&解凍。

URL="https://note.com/api/v2/attachments/download/513445af59f8bc1e9d9013ebd1f76fcb" && \
mkdir -p copilot_resample && cd copilot_resample && curl -L $URL | tar -xz

コンパイル方法

コンパイラの最適化はよくわかってないです
何かあればCopilotへ

nvcc -O3 -o resample resample.cpp resample.cu -Xcompiler="-O3 -march=native -fno-fast-math -funroll-loops -fopenmp" -lcudart -lsndfile -lfftw3

実行

実行方法が少し変わりました
非整数倍に対応してます
現時点ではwav(32bit)とflac(24bit)の対応

./resample "input.wav" "output.wav" -s 384000 -b 32 -sl 1024 -kb 10 -lp 21000 -th 16
入力ファイル            libsndfileが対応している拡張子
出力ファイル            wav(16~32bit)、flac(16~24bit)

-s  出力サンプリングレート    (数値指定で周波数、xと数値で係数、省略時はx2)
-b  出力ビット深度        (省略時、固定小数点の最大値)
-sl sinc_len              (省略時、自動計算)
-kb カイザー窓のβ         (省略時、自動計算)
-lp ローパスフィルターの追加処理 (省略時、無効)
-th マルチスレッド数       (省略時、自動計算)

最小構成はこれ

./resample "input.wav" "output.wav"

-lpモードについて

Copilot君と議論したところ、無限長のSinc関数は理想的なローパスフィルターとして動作するらしいです。
今回、膨大なSin_Lenを利用するように設計しているため、「半」理想的なローパスフィルターにぐらいはなりそうです。
そこで"-lp"モードを追加し、Sinc関数によるリサンプリングを2回実行するようにしました。
1回目をローパスフィルターとして使い、2回目をアップサンプリングに使います。
-lpで指定する値はローパスフィルターの遮断周波数ですが、実体としては1回目のリサンプリング時のナイキスト周波数です。
そのため、-lpで指定した2倍の周波数にリサンプリングしているだけです。

./resample "input.wav" "output.wav" -s 44100 -lp 22050
./resample "input.wav" "output.wav" -s 44100

この2つの指定は等価です。
2回同じサンプリングレートに変換するのは無駄なので、サンプリングレートが一致する場合は自動的に1回の処理になるようにしています。
アップサンプリングの一例ですが、以下のようにすると44100Hzから直接384000Hzにアップサンプリングするより高速になり、さらに「半」理想的なローパスフィルターまで適用されます。

./resample "input.wav" "output.wav" -s 384000 -lp 21000

21000Hzでローパスをしながら(実態は42000Hzへのリサンプリング)、44100Hz系から48000Hz系への変換が高速になります。(5分の音源で1日→1時間くらい違う)
サンプリング周波数の最大公約数を大きくするのが高速化のコツです。

Copilot君の激押しっぷりを見ると、DACのアンチエイリアスフィルターよりは高性能かも?


メモ

低音しっかり低い
高音いろんな音色する
高調波ガビガビしない
音がピタッと止まってリンギングない
問題なさそう

初めてC++触ったけど、ポインタとかメモリ触れるのに幸せを感じた

グラボよりFPGAの方がいいの?って何度も気になった
うたい文句通り数十倍速いなら、整数倍ならリアルタイム再生が視野に入る
でもそれだとYouTubeが…

12時間のハイレゾフォーマットを保存するための簡単な方法が無かった

ラズパイ5ってやっぱり進化の方向間違えた感ある


コードの内容

resample.cu

#include <cuda_runtime.h>
#include <iostream>
#include <cassert>

__device__ long long atomicAddl(unsigned long long* address, unsigned long long val) {
    unsigned long long old = *address, assumed;

    do {
        assumed = old;
        old = atomicCAS(address, assumed, val + assumed);
    } while (assumed != old);

    return old;
}

__device__ bool safeCast_longLong(const double& value, long long* result) {
    if (value < 0.0 || value > static_cast<double>(LLONG_MAX)) {
        return true;
    } else {
        *result = static_cast<long long>(value);
        return false;
    }
}

__device__ bool check_mul_overflow(const long long& a, const long long& b, long long* result) {

    if (a == 0 || b == 0) {
        *result = 0LL;
        return false;
    } else if (a > ULLONG_MAX / b) {
        return true;
    } else {
        *result = a * b;
        return false;
    }
}

__global__ void resampleKernel(double* inputData, double* outputData, long long inputChannels, long long inputFrames, long long outputFrames, long long sincLen, long long cycleFactor_Input, long long cycleFactor_Output, double* sincTable, unsigned long long* progressCounter, int* stopFlag) {
    const long long t = static_cast<long long>(blockIdx.x * blockDim.x + threadIdx.x);

    if (t >= outputFrames || *stopFlag != 0) {
        return;
    }

    long long temp_calcValue;
    if( check_mul_overflow(t, cycleFactor_Input, &temp_calcValue) || safeCast_longLong(floor(static_cast<double>(temp_calcValue) / static_cast<double>(cycleFactor_Output)), &temp_calcValue) ) {
        if(*stopFlag == 0) {
            atomicAdd(stopFlag, 1);
            printf("\nError CUDA: 計算時にオーバーフローしました");
        }
        assert(false);
    }

    const long long inputPos = temp_calcValue;
    const long long i        = t % cycleFactor_Output;

    for (long long c = 0; c < inputChannels; ++c) {
        double sum        = 0.0;
        const long long N = 2 * sincLen + 1;

        for (long long k = - sincLen; k <= sincLen; ++k) {
            if( *stopFlag != 0 ) {
                return;
            }

            const long long targetPos = inputPos + k;

            if (targetPos >= 0 && targetPos < inputFrames) {
                sum += inputData[targetPos * inputChannels + c] * sincTable[i * N + k + sincLen];
            }
        }
        outputData[t * inputChannels + c] = sum;
    }

    const unsigned long long progressInterval = static_cast<unsigned long long>(8LL * cycleFactor_Output);

    if (t + progressInterval <= outputFrames) {
        if ( t % progressInterval == 0LL ) {
            atomicAddl(progressCounter, progressInterval);
        }
    } else {
        atomicAddl(progressCounter, 1ULL);
    }
}

extern "C" void callResampleKernel(double* inputData, double* outputData, long long inputChannels, long long inputFrames, long long outputFrames, long long sincLen, long long cycleFactor_Input, long long cycleFactor_Output, double* sincTable, unsigned long long* progressCounter, int* stopFlag, cudaStream_t kernelStream) {

    const auto safeCast_unsignedInt = [](const long long& value){
        if ( value < 0 ) {
            throw std::overflow_error("負の値を変換できません");
        } else if ( value > static_cast<long long>(std::numeric_limits<unsigned int>::max())) {
            throw std::overflow_error("計算時にオーバーフローが発生しました");
        }

        return static_cast<unsigned int>(value);
    };

    const long long    totalThreads       = outputFrames;
  //const long long    numSM              =           28;
  //const long long    numCoresPerSM      =          128;
    const long long    maxThreadsPerBlock =         1024;

    const long long    totalBlocks        = (totalThreads + (maxThreadsPerBlock - 1LL)) / maxThreadsPerBlock;

    const unsigned int numBlocks          = safeCast_unsignedInt(totalBlocks);
    const unsigned int numThreadsPerBlock = safeCast_unsignedInt(maxThreadsPerBlock);

    resampleKernel<<<numBlocks, numThreadsPerBlock, 0 , kernelStream>>>(inputData, outputData, inputChannels, inputFrames, outputFrames, sincLen, cycleFactor_Input, cycleFactor_Output, sincTable, progressCounter, stopFlag );
}

resample.cpp

#include <iostream>
#include <span>
#include <vector>
#include <cmath>
#include <limits>
#include <string>
#include <cstring>
#include <sndfile.h>
#include <thread>
#include <iomanip>
#include <chrono>
#include <algorithm>
#include <map>
#include <functional>
#include <omp.h>
#include <fftw3.h>
#include <cuda_runtime.h>
#include "resample_error.h"
#include "resample_memma.h"
#include "resample_safecalc.h"

extern "C" void callResampleKernel(double* inputData, double* outputData, long long inputChannels, long long inputFrames, long long outputFrames, long long sincLen, long long cycleFactor_Input, long long cycleFactor_Output, double* sincTable, unsigned long long* progressCounter, int* d_stopFlag, cudaStream_t kernelStream);

class KaiserWindow {
public:
    void calcParam(const long long& inputSamplerate, const unsigned int& outputBitDepth, long long& sincLen, double& beta) {
        KaiserWindow kw;
        MEMManager& memma = MEMManager::getInstance();

        const double       dynamicRange = 20.0 * std::log10(std::pow(2, outputBitDepth) - 1);
        const long long    temp_sincLen = inputSamplerate;
        const long long    windowSize   = safeCalc(2LL) * temp_sincLen + 1LL;
        const long long    fftScale     = 500; // 約0.1Hzの精度
        const long long    fftSize      = ( safeCalc(2LL) * temp_sincLen ) * fftScale;
        const size_t       fftComplexAllocateSize = safeCalc(sizeof(fftw_complex)) * safeCalc(fftSize).to_size_t();


        fftw_complex *in  = (fftw_complex*) fftw_malloc(fftComplexAllocateSize);
        if(!in) {
            throw CustomError(ErrorType::FFTInitialization, "メモリの確保に失敗しました。」");
        }
        memma.store(in);

        fftw_complex *out = (fftw_complex*) fftw_malloc(fftComplexAllocateSize);
        if(!out) {
            throw CustomError(ErrorType::FFTInitialization, "メモリの確保に失敗しました。」");
        }
        memma.store(out);

        fftw_plan p = fftw_plan_dft_1d(fftSize, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
        if(!p) {
            throw CustomError(ErrorType::FFTInitialization, "メモリの確保に失敗しました。」");
        }
        memma.store(p);

        memset(in,  0, fftComplexAllocateSize);
        memset(out, 0, fftComplexAllocateSize);

        std::vector<long double> KWTable        (windowSize);
        std::vector<double>      KWTableFFTSpec (fftSize);
        std::vector<double>      KWTableFFT_DBFS(fftSize);

        const bool calcBeta    =    beta == 0.0 ? true : false;
        const bool calcSincLen = sincLen == 0LL ? true : false;
        
        if(calcBeta || calcSincLen)
        {
            if(calcBeta) beta =7.0;
            bool calcDigit = false;

            for(double loop = calcBeta ? 100.0 : (beta + 1.0); beta <= loop; beta += calcDigit ? 0.1 : 1.0 )
            {

                KWTable = kw.precomputed(temp_sincLen, beta);

                #pragma omp parallel for
                for (long long i = 0; i < windowSize; ++i) {
                    in[i][0] = safeCalc(KWTable[i]).to_double();
                    in[i][1] = 0.0;
                }

                fftw_execute(p);

                #pragma omp parallel for
                for (long long i = 0; i < fftSize; ++i) {
                    KWTableFFTSpec[i] = sqrt(safeCalc(out[i][0]) * out[i][0] + safeCalc(out[i][1]) * out[i][1]);
                }

                KWTableFFT_DBFS = kw.convertToDBFS(KWTableFFTSpec, fftSize);

                const long long secondPeakIndex = kw.findSecondPeakIndex(KWTableFFT_DBFS, fftSize);

                if(calcBeta) {
                    if( secondPeakIndex == -1LL || -dynamicRange < KWTableFFT_DBFS[secondPeakIndex] ) {
                        if( beta < loop ) continue;
                        else {
                            throw CustomError(ErrorType::CalcError, "適切なカイザー窓のβを見つけられませんでした。");
                        }
                    }

                    if(!calcDigit) {
                        calcDigit = true;
                        --beta;
                        continue;
                    }
                } else {
                    if( secondPeakIndex == -1LL ) {
                        throw CustomError(ErrorType::CalcError, "適切なカイザー窓のβを見つけられませんでした。より大きなβの値を試してみてください。");
                    }
                }

                std::cout << "β = " << beta << "のサイドローブのピーク (dB): " << KWTableFFT_DBFS[secondPeakIndex] << std::endl;

                if(calcSincLen) {
                    const long long sincLenFactorIndex = kw.findSincLenFactorIndex(KWTableFFT_DBFS, secondPeakIndex);

                    if(sincLenFactorIndex == -1LL) {
                        throw CustomError(ErrorType::CalcError, "適切なSinc_Lenを見つけられませんでした");
                    }

                    const long long marginScale        = 2LL;
                    const long long mainlobeScalePerHz = (safeCalc(sincLenFactorIndex) + (fftScale - 1LL)) * marginScale / fftScale;
                    sincLen = mainlobeScalePerHz * temp_sincLen;
                    std::cout << "1Hz未満の分解能を持つSinc_Lenの値: " << sincLen << std::endl;
                }

                break;
            }
        }

        memma.release(p);
        memma.release(in);
        memma.release(out);
    }

    std::vector<long double> precomputed(const long long sincLen, const double beta){

        KaiserWindow kw;

        const long long windowSize = safeCalc(2LL) * sincLen + 1LL;

        std::vector<long double> KWTable(windowSize);

        const long double collectWindowGain = kw.calcValue(sincLen, windowSize, beta);

        #pragma omp parallel for
        for (long long k = - sincLen; k <= sincLen; ++k) {
            KWTable[k + sincLen] = kw.calcValue(k + sincLen, windowSize, beta) / collectWindowGain;
        }

        KWTable[sincLen] = 1.0L;

        return KWTable;
    }

private:
    long double calcValue(const long long n, const long long N, const double beta) const {
        const long double x = ( (safeCalc(2LL) * n).to_long_double() / (safeCalc(N) - 1LL).to_long_double() ) - 1.0L;
        return std::cyl_bessel_il(0, safeCalc(beta).to_long_double() * sqrtl(safeCalc(1.0L) - powl(x, 2.0L))) / std::cyl_bessel_il(0, safeCalc(beta).to_long_double());
    }

    std::vector<double> convertToDBFS(const std::vector<double>& fftData, const long long& fftSize) const {
        std::vector<double> dbFSData(fftSize);

        const double maxAmplitude = *std::max_element(fftData.begin(), fftData.end());

        #pragma omp parallel for
        for (long long i = 0; i < fftSize; ++i) {
            const double normalizedValue = safeCalc(fftData[i]) / maxAmplitude;
            dbFSData[i]                  = safeCalc(20.0) * std::log10(normalizedValue);
        }

        return dbFSData;
    }

    long long findSecondPeakIndex(const std::vector<double>& KWTableFFTSpec, const long long& fftSize) const {
        const std::vector<double>::const_iterator mainLobePeak = std::max_element(KWTableFFTSpec.begin(), KWTableFFTSpec.end());
        const long long mainLobePeakIndex = safeCalc(std::distance(KWTableFFTSpec.begin(), mainLobePeak)).to_long_long();
              bool foundDescent = true;

        for (long long i = mainLobePeakIndex + 1LL; i < fftSize - 1LL; ++i) {
            if (foundDescent) {
                if (KWTableFFTSpec[i] < KWTableFFTSpec[i + 1LL]) {
                    foundDescent = false;
                }
            } else {
                if (KWTableFFTSpec[i] > KWTableFFTSpec[i + 1LL]) {
                    return i;
                }
            }
        }

        return -1LL;
    }

    long long findSincLenFactorIndex(const std::vector<double>& KWTableFFTSpec, const long long& secondPeakIndex) const {
        const double secondLobePeakValue = KWTableFFTSpec[secondPeakIndex];

        for (long long i = secondPeakIndex - 1LL; i > 0LL; --i) {
            if (KWTableFFTSpec[i] > secondLobePeakValue) {
                return i + 1LL;
            }
        }

        return -1LL;
    }
};

class Sinc {

public:
    std::vector<std::vector<double>> precomputed(const long long sincLen, const long long& cycleFactor_Input, const long long& cycleFactor_Output, const std::vector<long double>& kaiserWindowTable) {

        const double oversamplingFactor = safeCalc(cycleFactor_Output).to_double() / safeCalc(cycleFactor_Input).to_double();

        std::vector<std::vector<double>> sincTable(cycleFactor_Output, std::vector<double>(safeCalc(2LL) * sincLen + 1LL, 0.0));

        const auto calcSampleGap = [&cycleFactor_Input, &cycleFactor_Output](const long long &i) -> long double {
            const long long   temp_calcValue    = safeCalc(i) * cycleFactor_Input;
            const long double outPosMappedInPos = safeCalc(temp_calcValue).to_long_double() / safeCalc(cycleFactor_Output).to_long_double();

            return floorl(outPosMappedInPos) - outPosMappedInPos;
        };

        #pragma omp parallel for
        for (long long i = 0LL; i < cycleFactor_Output; ++i) {
            const long double sampleGap = calcSampleGap(i);

            for (long long k = - sincLen; k <= sincLen; ++k) {
                const long double value = calcValue(safeCalc(sampleGap) + safeCalc(k).to_long_double()) * kaiserWindowTable[k + sincLen];
                sincTable[i][k + sincLen] = safeCalc(value).to_double();
            }
        }
        return sincTable;
    }

private:
    long double calcValue(const long double x) {
        if ( x == 0.0 ) {
            return 1.0L;
        } else {
            const long double y = safeCalc(M_PIl) * x;
            return sin(y) / y;
        }
    }

};

void resampleProgress(unsigned long long* ull_progressCounter, const long long& outputFrames, cudaStream_t& progressStream) {

          unsigned long long progressCounter     = 0;
    const double             startTime           = omp_get_wtime();
    const double             double_outputFrames = static_cast<double>(outputFrames);
    const unsigned long long ull_outputFrames    = static_cast<unsigned long long>(outputFrames);

    auto timeUnitOptmz = [](double t) -> std::string {
        const double SECONDS_IN_A_DAY    = 86400.0;
        const double SECONDS_IN_AN_HOUR  = 3600.0;
        const double SECONDS_IN_A_MINUTE = 60.0;

        std::string timeUnit;

        if (t >= SECONDS_IN_A_DAY   ) {
            timeUnit       = "日";
            t /= SECONDS_IN_A_DAY;
        } else if (t >= SECONDS_IN_AN_HOUR) {
            timeUnit       = "時間";
            t /= SECONDS_IN_AN_HOUR;
        } else if (t >= SECONDS_IN_A_MINUTE) {
            timeUnit       = "分";
            t /= SECONDS_IN_A_MINUTE;
        } else {
            timeUnit = "秒";
        }

        std::ostringstream oss;
        oss << std::fixed << std::setprecision(2) << t;

        return oss.str() + " " + timeUnit;
    };

    while (progressCounter < ull_outputFrames) {
        cudaMemcpyAsync(&progressCounter, ull_progressCounter, sizeof(unsigned long long), cudaMemcpyDeviceToHost, progressStream);

        cudaError_t err = cudaGetLastError();
        if (err != cudaSuccess) {
            return;
        }
        cudaStreamSynchronize(progressStream);

        const double elapsedTime = omp_get_wtime() - startTime;
        if (progressCounter >= outputFrames) break;
        if (progressCounter >  0ULL) {
            const double progress           = 100.0       * (static_cast<double>(progressCounter) / double_outputFrames);
            const double estimatedTotalTime = elapsedTime * (double_outputFrames / static_cast<double>(progressCounter));
            const double remainingTime      = estimatedTotalTime - elapsedTime;

            std::cout << "\r進行中: " << std::fixed << std::setprecision(2) << progress << " %"
                      << " 残り時間: " << timeUnitOptmz(remainingTime) << std::flush;
        }
        std::this_thread::sleep_for(std::chrono::seconds(1));
    }

    const double totalTime = omp_get_wtime() - startTime;
    std::cout << "\rリサンプリング進行中: 100.00 % 経過時間: " << timeUnitOptmz(totalTime) << std::endl << std::flush;

}

void resampleWithCUDA(double* inputData, double* outputData, const long long& inputChannels, const long long& inputFrames, const long long& outputFrames, const long long& sincLen, const long long& cycleFactor_Input, const long long& cycleFactor_Output, const std::vector<std::vector<double>>& sincTable) {
    MEMManager& memma = MEMManager::getInstance();

    double*             d_inputData;
    double*             d_outputData;
    double*             d_sincTable;
    unsigned long long* ull_progressCounter;
    int*                d_stopFlag;
    const size_t inputDataAllocSize  = safeCalc(inputChannels).to_size_t() * safeCalc(inputFrames ).to_size_t() * sizeof(double);
    const size_t outputDataAllocSize = safeCalc(inputChannels).to_size_t() * safeCalc(outputFrames).to_size_t() * sizeof(double);
    size_t cudaAllocMaxSize          = safeCalc(std::numeric_limits<unsigned int>::max()).to_size_t();

    std::vector<double> flattenedSincTable;
    size_t flattenedSincTableAllocSize;


    for (const auto& row : sincTable) {
        flattenedSincTable.insert(flattenedSincTable.end(), row.begin(), row.end());
    }

    flattenedSincTableAllocSize = safeCalc(flattenedSincTable.size()) * sizeof(double);

    const auto cudaSafeMalloc = [&memma](auto*& pointer, size_t allocSize) {
        const cudaError_t err = cudaMalloc((void**)&pointer,  allocSize);
        if (err != cudaSuccess) {
            throw CustomError(ErrorType::CUDAInitialization, "メモリの初期化に失敗");
        }

        memma.store(pointer);

        cudaMemset(pointer, 0, allocSize);
        if (err != cudaSuccess) {
            throw CustomError(ErrorType::CUDAInitialization, "メモリの初期化に失敗");
        }
    };

    const auto cudaSafeMallocManaged = [&memma](auto*& pointer, size_t allocSize) {
        const cudaError_t err = cudaMallocManaged((void**)&pointer,  allocSize);
        if (err != cudaSuccess) {
            throw CustomError(ErrorType::CUDAInitialization, "メモリの初期化に失敗");
        }

        memma.store(pointer);

        cudaMemset(pointer, 0, allocSize);
        if (err != cudaSuccess) {
            throw CustomError(ErrorType::CUDAInitialization, "メモリの初期化に失敗");
        }
    };


    const auto cudaSafeMemcpy = [](auto*& pointer, auto* data, const auto& allocSize, const cudaMemcpyKind& kind) {
        const cudaError_t err = cudaMemcpy(pointer, data, allocSize, kind);
        if (err != cudaSuccess) {
            throw CustomError(ErrorType::CUDAInitialization, "メモリの初期化に失敗");
        }
    };

    const auto cudaSafeStreamCreate = [&memma](cudaStream_t& stream) {
        const cudaError_t err = cudaStreamCreate(&stream);
        if (err != cudaSuccess) {
            throw CustomError(ErrorType::CudaError, "streamの作成に失敗: " + *cudaGetErrorString(err));
        }
        memma.store(stream);
    };

    cudaStream_t kernelStream;
    cudaStream_t progressStream;

    cudaSafeStreamCreate(kernelStream);
    cudaSafeStreamCreate(progressStream);

    const auto stopCUDAKernel = [&d_stopFlag, &memma, &progressStream]() {
            int stop = 1;
            cudaMemcpyAsync(d_stopFlag, &stop, sizeof(int), cudaMemcpyHostToDevice, progressStream);
            cudaStreamSynchronize(progressStream);
    };
    memma.setExitFunction(stopCUDAKernel);

    if (inputDataAllocSize > cudaAllocMaxSize) {
        cudaSafeMallocManaged(d_inputData,  inputDataAllocSize);
    } else {
        cudaSafeMalloc(d_inputData,  inputDataAllocSize);
    }


    if (outputDataAllocSize > cudaAllocMaxSize) {
        cudaSafeMallocManaged(d_outputData, outputDataAllocSize);
    } else {
        cudaSafeMalloc(d_outputData, outputDataAllocSize);
    }

    if (flattenedSincTableAllocSize > cudaAllocMaxSize) {
        cudaSafeMallocManaged(d_sincTable, flattenedSincTableAllocSize);
    } else {
        cudaSafeMalloc(d_sincTable, flattenedSincTableAllocSize);
    }

    cudaSafeMalloc(ull_progressCounter, sizeof(unsigned long long));
    cudaSafeMalloc(d_stopFlag, sizeof(int));

    cudaSafeMemcpy(d_inputData, inputData, inputDataAllocSize, cudaMemcpyHostToDevice);
    cudaSafeMemcpy(d_sincTable, flattenedSincTable.data(), flattenedSincTableAllocSize, cudaMemcpyHostToDevice);

    unsigned long long initialCounter = 0;
    cudaSafeMemcpy(ull_progressCounter, &initialCounter, sizeof(unsigned long long), cudaMemcpyHostToDevice);

    int initialFlag = 0;
    cudaSafeMemcpy(d_stopFlag, &initialFlag, sizeof(int), cudaMemcpyHostToDevice);

    cudaDeviceSynchronize();

    callResampleKernel(d_inputData, d_outputData, inputChannels, inputFrames, outputFrames, sincLen, cycleFactor_Input, cycleFactor_Output, d_sincTable, ull_progressCounter, d_stopFlag, kernelStream);
    resampleProgress(ull_progressCounter, outputFrames, progressStream);

    cudaStreamSynchronize(kernelStream);
    cudaSafeMemcpy(outputData, d_outputData, safeCalc(inputChannels).to_size_t() * safeCalc(outputFrames).to_size_t() * sizeof(double), cudaMemcpyDeviceToHost);

    const cudaError_t err = cudaGetLastError();
    if (err != cudaSuccess) {
        throw CustomError(ErrorType::CudaError, std::string(cudaGetErrorString(err)));
    }

    memma.releaseAll();
}

int getLibsndfileFormat (const std::string& audioFormat, unsigned int& outputBitDepth) {

    int format;

    if (audioFormat == "WAV") {
        if (outputBitDepth == 0) {
            outputBitDepth = 32;
        }
        format = outputBitDepth == 32 ? SF_FORMAT_WAV  | SF_FORMAT_PCM_32 :
                 outputBitDepth == 24 ? SF_FORMAT_WAV  | SF_FORMAT_PCM_24 :
                 outputBitDepth == 16 ? SF_FORMAT_WAV  | SF_FORMAT_PCM_16 : -1;
    } else if (audioFormat == "FLAC") {
        if (outputBitDepth == 0) {
            outputBitDepth = 24;
        }
        format = outputBitDepth == 24 ? SF_FORMAT_FLAC | SF_FORMAT_PCM_24 :
                 outputBitDepth == 16 ? SF_FORMAT_FLAC | SF_FORMAT_PCM_16 : -1;
    } else {
        format = -1;
    }

    if (format != -1) {
        return format;
    } else {
        throw CustomError(ErrorType::InvalidParameter, "出力フォーマットが非対応です");
    }
}

void resampleAudio(const std::string& inputFile, const std::string& outputFile, long long & outputSamplerate, unsigned int& outputBitDepth, const long double& oversamplingFactor, long long& sincLen, double& beta, const long long& lowPass, const unsigned int& numThreads) {


    std::cout << "データ読み込み中..." << std::endl;


    SF_INFO sfinfo;
    SNDFILE *infile = sf_open(inputFile.c_str(), SFM_READ, &sfinfo);
    if (!infile) {
        throw CustomError(ErrorType::FileNotOpen, "入力ファイルを開けませんでした(" + inputFile + ") " + sf_strerror(NULL));
    }

    const auto nextMultiple = [](long long x, long long scale) -> long long {
        return ((x + (scale -1))/scale) * scale;
    };

    const long long     inputSamplerate = safeCalc(sfinfo.samplerate).to_long_long();
    const long long     inputFrames     = nextMultiple(safeCalc(sfinfo.frames).to_long_long(), inputSamplerate); // オリジナルの長さを破壊してx.0秒にゼロパディングする動作に統一
    const long long     inputChannels   = safeCalc(sfinfo.channels).to_long_long();
    std::vector<double> inputData(inputFrames * inputChannels, 0.0);
    sf_read_double(infile, inputData.data(), inputData.size());
    sf_close(infile);


    std::cout << "書き出し設定の事前確認中..." << std::endl;


    SF_INFO outsfinfo;
    if (oversamplingFactor != 0.0L) {
        outputSamplerate = safeCalc(floorl(static_cast<long double>(inputSamplerate) * oversamplingFactor)).to_long_long();
    }

    std::string audioFormat    = "";
    std::string outputBaseName = "";
    std::string outputExt      = "";
    [&outputFile, &audioFormat, &outputBaseName, &outputExt]() {
        std::string ext4 = outputFile.substr(outputFile.length() - 4);
        std::string ext5 = outputFile.substr(outputFile.length() - 5);

        if (ext4 == ".wav") {
            audioFormat    = "WAV";
            outputBaseName = outputFile.substr(0, outputFile.length() - 4);
            outputExt      = ext4.substr(1);       
        } else if (ext5 == ".flac") {
            audioFormat    = "FLAC";
            outputBaseName = outputFile.substr(0, outputFile.length() - 5);
            outputExt      = ext5.substr(1);     
        }
    }();

    outsfinfo.samplerate = safeCalc(outputSamplerate).to_int();
    outsfinfo.channels   = safeCalc(inputChannels).to_int();
    outsfinfo.format     = getLibsndfileFormat(audioFormat, outputBitDepth);


    // オリジナルの長さを破壊してx.0秒にゼロパディングする動作に統一
    const long long outputFrames = ((safeCalc(inputFrames) + inputSamplerate - 1LL) / inputSamplerate) * outputSamplerate;
    const long long outputSize   = safeCalc(outputFrames) * inputChannels;


    long long cycleFactor_Input;
    long long cycleFactor_Output;

    const auto calcCycleFactor = [](const long long& inputSamplerate, const long long& outputSamplerate, long long& cycleFactor_Input, long long& cycleFactor_Output) {
        const long long samplerate_gcd = [](long long x, long long y) -> long long {
            while (y != 0) {
                long long temp = y;
                y = x % y;
                x = temp;
            }
            return x;
        }(inputSamplerate, outputSamplerate);
        cycleFactor_Input  = inputSamplerate  / samplerate_gcd;
        cycleFactor_Output = outputSamplerate / samplerate_gcd;
    };

    std::vector<double> outputData(outputSize, 0.0);


    KaiserWindow kw;
    Sinc sinc;
    std::vector<long double> kaiserWindowTable;
    std::vector<std::vector<double>> sincTable;

    if (lowPass > 0LL && lowPass != outputSamplerate / 2LL) {
        std::vector<double> lowPassedData(outputSize, 0.0);

        double    lowPass_beta         = beta;
        long long lowPass_sincLen      = sincLen;
        long long lowPass_samplerate   = safeCalc(lowPass) * 2LL;
        long long lowPass_frames = ((safeCalc(inputFrames) + inputSamplerate - 1LL) / inputSamplerate) * lowPass_samplerate;
        long long lowPass_cycleFactor_Input;
        long long lowPass_cycleFactor_Output;

        calcCycleFactor(inputSamplerate, lowPass_samplerate, lowPass_cycleFactor_Input, lowPass_cycleFactor_Output);


        std::cout << "ローパス中..." << std::endl;


        kw.calcParam(static_cast<long long>(inputSamplerate), outputBitDepth, lowPass_sincLen, lowPass_beta);
        kaiserWindowTable = kw.precomputed(lowPass_sincLen, lowPass_beta);
        sincTable         = sinc.precomputed(lowPass_sincLen, lowPass_cycleFactor_Input, lowPass_cycleFactor_Output, kaiserWindowTable);

        resampleWithCUDA(inputData.data(), lowPassedData.data(), inputChannels, inputFrames, lowPass_frames, lowPass_sincLen, lowPass_cycleFactor_Input, lowPass_cycleFactor_Output, sincTable);
        
        calcCycleFactor(lowPass_samplerate, outputSamplerate, cycleFactor_Input, cycleFactor_Output);


        std::cout << "リサンプリング中..." << std::endl;


        kw.calcParam(lowPass_samplerate, outputBitDepth, sincLen, beta);
        kaiserWindowTable = kw.precomputed(sincLen, beta);
        sincTable         = sinc.precomputed(sincLen, cycleFactor_Input, cycleFactor_Output, kaiserWindowTable);

        resampleWithCUDA(lowPassedData.data(), outputData.data(), inputChannels, lowPass_frames, outputFrames, sincLen, cycleFactor_Input, cycleFactor_Output, sincTable);
    } else {
        calcCycleFactor(inputSamplerate, outputSamplerate, cycleFactor_Input, cycleFactor_Output);


        std::cout << "リサンプリング中..." << std::endl;


        kw.calcParam(inputSamplerate, outputBitDepth, sincLen, beta);
        kaiserWindowTable = kw.precomputed(sincLen, beta);
        sincTable         = sinc.precomputed(sincLen, cycleFactor_Input, cycleFactor_Output, kaiserWindowTable);

        resampleWithCUDA(inputData.data(), outputData.data(), inputChannels, inputFrames, outputFrames, sincLen, cycleFactor_Input, cycleFactor_Output, sincTable);
    }


    std::cout << "正規化中..." << std::endl;


    double maxAmplitude = 0.0;

    #pragma omp parallel for
    for (double sample : outputData) {
        const double absSample = std::abs(sample);
        #pragma omp critical
        {
            if (absSample > maxAmplitude) {
                maxAmplitude = absSample;
            }
        }
    }

    if (maxAmplitude > 0.0) {
        #pragma omp parallel for
        for (double& sample : outputData) {
            sample /= maxAmplitude;
        }
    }


    std::cout << "データ書き込み中..." << std::endl;


    auto formatNumberWithDigits = [](unsigned long long number, unsigned long long digits) -> std::string {
        std::ostringstream oss;
        oss << std::setw(digits) << std::setfill('0') << number;
        return oss.str();
    };

    const unsigned long long maxChunks   = []() -> unsigned long long {
        const unsigned long long maxFileSize = safeCalc(std::numeric_limits<std::uint32_t>::max()).to_unsigned_long_long();
        const unsigned long long headerSize  = 44ULL;
        return (maxFileSize - headerSize ) / safeCalc(sizeof(uint32_t)).to_unsigned_long_long();
    }();

          unsigned long long fileIndex          = 0ULL;
    const unsigned long long outputDataSize     = safeCalc(outputData.size()).to_unsigned_long_long();
    const unsigned long long maxFileCount       = (outputDataSize + maxChunks - 1ULL) / maxChunks;
    const unsigned long long maxFileCountDigits = safeCalc(std::to_string(maxFileCount).length()).to_unsigned_long_long();

    MEMManager& memma = MEMManager::getInstance();

    for (unsigned long long startChunk = 0ULL; startChunk < outputDataSize; startChunk += maxChunks) {
        const unsigned long long endChunk = safeCalc(startChunk) + maxChunks;

        const std::vector<double> subset(outputData.begin() + startChunk, (endChunk > outputDataSize) ? outputData.end() : (outputData.begin() + endChunk));

        std::string outputFileName = outputBaseName;
        if( maxFileCount > 1ULL ) {
            outputFileName += "_" + formatNumberWithDigits(++fileIndex, maxFileCountDigits);
        }
        outputFileName += "." + outputExt;

        SNDFILE *outfile = sf_open(outputFileName.c_str(), SFM_WRITE, &outsfinfo);
        if (!outfile) {
            throw CustomError(ErrorType::FileNotOpen, "出力ファイルを開けませんでした(" + outputFileName + ") " + sf_strerror(NULL));
        }
        memma.store(outfile);

        sf_write_double(outfile, subset.data(), subset.size());
        memma.release(outfile);
    }


    std::cout << "リサンプリング完了!" << std::endl;
}

class OptionHandler {
public:
    static void processS(const std::string &arg, long double& oversamplingFactor, long long& outputSamplerate) {
        if (arg.substr(0, 1) == "x") {
            oversamplingFactor = safeCalc(std::stod(arg.substr(1))).to_long_double();
        } else {
            outputSamplerate = std::stoll(arg, nullptr, 10);
            if (outputSamplerate < 0LL) {
                std::cerr << "出力サンプリング周波数に指定できるのは正値のみです" << std::endl;
                exit(1);
            }
        }
    }

    static void processB(const std::string &arg, unsigned int& outputBitDepth) {
        outputBitDepth = static_cast<unsigned int>(std::atoi(arg.c_str()));
    }

    static void processSL(const std::string &arg, long long& sincLen) {
        sincLen = std::stoll(arg, nullptr, 10);
        if (sincLen < 0LL) {
            std::cerr << "Sinc_Lenに指定できるのは正値のみです" << std::endl;
            exit(1);
        }
    }

    static void processKB(const std::string &arg, double& beta) {
        beta = std::stod(arg);
    }

    static void processLP(const std::string &arg, long long& lowPass) {
        lowPass = std::stoll(arg, nullptr, 10);
        if (lowPass < 0LL) {
            std::cerr << "ローパスフィルターに指定できるのは正値のみです" << std::endl;
            exit(1);
        }
    }

    static void processTH(const std::string &arg, unsigned int& numThreads) {
        numThreads = static_cast<unsigned int>(std::atoi(arg.c_str()));
    }
};

int main(int argc, char *argv[]) {

    if (argc < 2) {
        std::cerr << "使用方法: " << argv[0] << " <出力ファイルパス> <入力ファイルパス> -s <出力サンプリング周波数:数値 数値の頭にxをつけると係数> -b <出力ビット深度:数値> -sl <Sinc_Len:数値> -kb <カイザー窓 β:小数値> -lp <ローパスフィルター周波数:数値> -th <並列処理数:数値>" << std::endl;
        return 1;
    }

    const std::string inputFile  = argv[1];
    const std::string outputFile = argv[2];

    double       beta               =  0.0;
    long double  oversamplingFactor = 0.0L;
    long long    outputSamplerate   =  0LL;
    long long    sincLen            =  0LL;
    unsigned int outputBitDepth     =   0U;
    long long    lowPass            =  0LL;
    unsigned int numThreads         =   0U;

    std::map<std::string, std::function<void(const std::string&)>> handlerMap = {
        {"-s" , [&](const std::string &arg) { OptionHandler::processS (arg, oversamplingFactor, outputSamplerate); }},
        {"-b" , [&](const std::string &arg) { OptionHandler::processB (arg, outputBitDepth); }},
        {"-sl", [&](const std::string &arg) { OptionHandler::processSL(arg, sincLen); }},
        {"-kb", [&](const std::string &arg) { OptionHandler::processKB(arg, beta); }},
        {"-lp", [&](const std::string &arg) { OptionHandler::processLP(arg, lowPass); }},
        {"-th", [&](const std::string &arg) { OptionHandler::processTH(arg, numThreads); }}
    };

    for (int i = 3; i < argc; i++) {
        auto it = handlerMap.find(argv[i]);
        if (it != handlerMap.end()) {
            if (i + 1 < argc && handlerMap.find(argv[i + 1]) == handlerMap.end()) {
                it->second(argv[++i]);
            } else {
                std::cerr << argv[i] << " にはパラメーターが必要です" << std::endl;
                exit(1);
            }
        }
    }

    try {
        if (numThreads != 0U) {
            omp_set_num_threads(numThreads);
        }
        resampleAudio(inputFile, outputFile, outputSamplerate, outputBitDepth, oversamplingFactor, sincLen, beta, lowPass, numThreads);
    } catch (const CustomError& e) {
        MEMManager& memma = MEMManager::getInstance();
        memma.releaseAll();

        std::cerr << e.what() << std::endl;
        return 1;
    }

    return 0;
}

resample_error.h

#ifndef RESAMPLE_ERROR_H
#define RESAMPLE_ERROR_H

#include <exception>
#include <string>

enum class ErrorType {
    InvalidParameter,
    FileNotOpen,
    FFTInitialization,
    CUDAInitialization,
    CudaError,
    CalcError,
    Unknown
};

class CustomError : public std::exception {
public:
    CustomError(ErrorType type, const std::string& details)
        : type_(type), message_(generateMessage(type, details)) {}

    const char* what() const noexcept override {
        return message_.c_str();
    }

private:
    ErrorType type_;
    std::string message_;

    std::string generateMessage(ErrorType type, const std::string& details) const {
        std::string prefix;
        switch (type) {
            case ErrorType::InvalidParameter:
                prefix = "Parameter Error: ";
                break;
            case ErrorType::FileNotOpen:
                prefix = "FileNotOpen Error: ";
                break;
            case ErrorType::FFTInitialization:
                prefix = "FFTInitialization Error: ";
                break;
            case ErrorType::CUDAInitialization:
                prefix = "CUDAInitialization Error: ";
                break;
            case ErrorType::CudaError:
                prefix = "CUDA Error: ";
                break;
            case ErrorType::CalcError:
                prefix = "Calc Error: ";
                break;
            default:
                prefix = "Unknown Error: ";
                break;
        }
        return prefix + details;
    }
};

#endif

resample_memma.h

#ifndef RESAMPLE_MEMMA_H
#define RESAMPLE_MEMMA_H

#include <csignal>
#include <fftw3.h>
#include <cuda_runtime.h>
#include <sndfile.h>
#include <mutex>

class MEMManager {
public:
    static MEMManager& getInstance() {
        static MEMManager instance;
        return instance;
    }

    void store(int* ptr) {
        std::lock_guard<std::mutex> lock(mtx);
        memory_pointers_cudaInt.push_back(ptr);
    }

    void store(double* ptr) {
        std::lock_guard<std::mutex> lock(mtx);
        memory_pointers_cudaDouble.push_back(ptr);
    }

    void store(unsigned int* ptr) {
        std::lock_guard<std::mutex> lock(mtx);
        memory_pointers_cudaUnsignedInt.push_back(ptr);
    }

    void store(long long* ptr) {
        std::lock_guard<std::mutex> lock(mtx);
        memory_pointers_cudaLongLong.push_back(ptr);
    }

    void store(unsigned long long* ptr) {
        std::lock_guard<std::mutex> lock(mtx);
        memory_pointers_cudaUnsignedLongLong.push_back(ptr);
    }

    void store(cudaStream_t ptr) {
        std::lock_guard<std::mutex> lock(mtx);
        memory_pointers_cudaStream.push_back(ptr);
    }

    void store(fftw_plan ptr) {
        std::lock_guard<std::mutex> lock(mtx);
        memory_pointers_plan.push_back(ptr);
    }

    void store(fftw_complex* ptr) {
        std::lock_guard<std::mutex> lock(mtx);
        memory_pointers_complex.push_back(ptr);
    }

    void store(SNDFILE* ptr) {

        std::lock_guard<std::mutex> lock(mtx);
        memory_pointers_sndfile.push_back(ptr);
    }

    template<typename T>
    void release(T* ptr, std::vector<T*>& container) {
        std::lock_guard<std::mutex> lock(mtx);
        auto it = std::find(container.begin(), container.end(), ptr);
        if (it != container.end()) {
            const cudaError_t err = cudaFree(*it);
            if (err != cudaSuccess) {
                std::cerr << "Failed to CUDA Free: " << cudaGetErrorString(err) << std::endl;
            }
            container.erase(it);
        }
    }

    void release(int* ptr) {
        release(ptr, memory_pointers_cudaInt);
    }

    void release(double* ptr) {
        release(ptr, memory_pointers_cudaDouble);
    }

    void release(unsigned int* ptr) {
        release(ptr, memory_pointers_cudaUnsignedInt);
    }

    void release(long long* ptr) {
        release(ptr, memory_pointers_cudaLongLong);
    }

    void release(unsigned long long* ptr) {
        release(ptr, memory_pointers_cudaUnsignedLongLong);
    }

    void release(cudaStream_t ptr) {
        std::lock_guard<std::mutex> lock(mtx);
        auto it = std::find(memory_pointers_cudaStream.begin(), memory_pointers_cudaStream.end(), ptr);
        if (it != memory_pointers_cudaStream.end()) {
            const cudaError_t err = cudaStreamDestroy(*it);
            if (err != cudaSuccess) {
                std::cerr << "Failed to Destroy Stream: " << cudaGetErrorString(err) << std::endl;
            }
            memory_pointers_cudaStream.erase(it);
        }
    }

    void release(fftw_plan ptr) {
        std::lock_guard<std::mutex> lock(mtx);
        auto it = std::find(memory_pointers_plan.begin(), memory_pointers_plan.end(), ptr);
        if (it != memory_pointers_plan.end()) {
            fftw_destroy_plan(*it);
            memory_pointers_plan.erase(it);
        }
    }

    void release(fftw_complex* ptr) {
        std::lock_guard<std::mutex> lock(mtx);
        auto it = std::find(memory_pointers_complex.begin(), memory_pointers_complex.end(), ptr);
        if (it != memory_pointers_complex.end()) {
            fftw_free(*it);
            memory_pointers_complex.erase(it);
        }
    }

    void release(SNDFILE* ptr) {
        std::lock_guard<std::mutex> lock(mtx);
        auto it = std::find(memory_pointers_sndfile.begin(), memory_pointers_sndfile.end(), ptr);
        if (it != memory_pointers_sndfile.end()) {
            sf_close(*it);
            memory_pointers_sndfile.erase(it);
        }
    }



    void releaseAll() {
        runExitFunction();

        for (auto ptr : memory_pointers_cudaStream) {
            cudaStreamDestroy(ptr);
        }
        memory_pointers_cudaStream.clear();

        for (auto ptr : memory_pointers_cudaInt) {
            cudaFree(ptr);
        }
        memory_pointers_cudaInt.clear();

        for (auto ptr : memory_pointers_cudaDouble) {
            cudaFree(ptr);
        }
        memory_pointers_cudaDouble.clear();

        for (auto ptr : memory_pointers_cudaUnsignedInt) {
            cudaFree(ptr);
        }
        memory_pointers_cudaUnsignedInt.clear();

        for (auto ptr : memory_pointers_cudaLongLong) {
            cudaFree(ptr);
        }
        memory_pointers_cudaLongLong.clear();

        for (auto ptr : memory_pointers_cudaUnsignedLongLong) {
            cudaFree(ptr);
        }
        memory_pointers_cudaUnsignedLongLong.clear();

        for (auto ptr : memory_pointers_plan) {
            fftw_destroy_plan(ptr);
        }
        memory_pointers_plan.clear();

        for (auto ptr : memory_pointers_complex) {
            fftw_free(ptr);
        }
        memory_pointers_complex.clear();

        for (auto ptr : memory_pointers_sndfile) {
            sf_close(ptr);
        }
        memory_pointers_sndfile.clear();
    }

    void setExitFunction(std::function<void()> func) {
        std::lock_guard<std::mutex> lock(mtx);
        exitFunction = func;
    }

    void runExitFunction() {
        if (exitFunction)
        {
            exitFunction();
        }
    }

    std::mutex& getMutex() {
        return mtx;
    }

private:
    std::mutex mtx;
    bool exitCalled = false;

    MEMManager() {
        signal(SIGINT, signalHandler);
    }

    ~MEMManager() {
        if (!exitCalled) {
            {
                std::lock_guard<std::mutex> lock(mtx);
                exitCalled = true;
            }
            releaseAll();
        }
    }

    MEMManager(const MEMManager&) = delete;
    MEMManager& operator=(const MEMManager&) = delete;

    static std::vector<int*>                memory_pointers_cudaInt;
    static std::vector<double*>             memory_pointers_cudaDouble;
    static std::vector<unsigned int*>       memory_pointers_cudaUnsignedInt;
    static std::vector<long long*>          memory_pointers_cudaLongLong;
    static std::vector<unsigned long long*> memory_pointers_cudaUnsignedLongLong;
    static std::vector<cudaStream_t>        memory_pointers_cudaStream;
    static std::vector<fftw_plan>           memory_pointers_plan;
    static std::vector<fftw_complex*>       memory_pointers_complex;
    static std::vector<SNDFILE*>            memory_pointers_sndfile;
    std::function<void()>                   exitFunction = nullptr;

    static void signalHandler(int signum) {
        MEMManager& memma = getInstance();
        if (!memma.exitCalled) {
            {
                std::lock_guard<std::mutex> lock(memma.mtx);
                memma.exitCalled = true;
            }
            std::cout << "\n処理が中断されました。" << std::endl;
            memma.releaseAll();
            exit(signum);
        }
    }
};


std::vector<int*>                MEMManager::memory_pointers_cudaInt;
std::vector<double*>             MEMManager::memory_pointers_cudaDouble;
std::vector<unsigned int*>       MEMManager::memory_pointers_cudaUnsignedInt;
std::vector<long long*>          MEMManager::memory_pointers_cudaLongLong;
std::vector<unsigned long long*> MEMManager::memory_pointers_cudaUnsignedLongLong;
std::vector<cudaStream_t>        MEMManager::memory_pointers_cudaStream;
std::vector<fftw_plan>           MEMManager::memory_pointers_plan;
std::vector<fftw_complex*>       MEMManager::memory_pointers_complex;
std::vector<SNDFILE*>            MEMManager::memory_pointers_sndfile;

#endif

resample_safecalc.h

#include <type_traits>
#include <limits>
#include <stdexcept>
#include <cmath>
#include <sndfile.h>

template <typename T>
class calcImpl {
public:
    T value;

    calcImpl(T v) : value(v) {}

    calcImpl operator+(T v) const {
        T result;
        if (is_integer()) {
            if (value > std::numeric_limits<T>::max() - v) {
                throw std::overflow_error("加算時にオーバーフローが発生しました");
            }
            result = value + v;
        } else {
            result = value + v;
            if (!std::isinf(value) && !std::isinf(v) && std::isinf(result)) {
                throw std::overflow_error("加算時にオーバーフローが発生しました");
            }
        }
        return calcImpl(result);
    }

    calcImpl operator-(T v) const {
        T result;
        if (is_integer()) {
            if (value < v + std::numeric_limits<T>::lowest()) {
                throw std::overflow_error("減算時にオーバーフローが発生しました");
            }
            result = value - v;
        } else {
            result = value - v;
            if (!std::isinf(value) && !std::isinf(v) && std::isinf(result)) {
                throw std::overflow_error("減算時にオーバーフローが発生しました");
            }
        }
        return calcImpl(result);
    }

    calcImpl operator*(T v) const {
        T result;
        if (is_integer()) {
            if (value != 0 && (v > std::numeric_limits<T>::max() / value || v < std::numeric_limits<T>::lowest() / value)) {
                throw std::overflow_error("乗算時にオーバーフローが発生しました");
            }
            result = value * v;
        } else {
            result = value * v;
            if (!std::isinf(value) && !std::isinf(v) && std::isinf(result)) {
                throw std::overflow_error("乗算時にオーバーフローが発生しました");
            }
        }
        return calcImpl(result);
    }

    calcImpl operator/(T v) const {
        if (v == 0) {
            throw std::invalid_argument("ゼロ除算は許可されていません");
        }
        return calcImpl(value / v);
    }

    template <typename U>
    calcImpl<U> to_type() const {
        check_overflow<U>(value);
        return calcImpl<U>(static_cast<U>(value));
    }

    calcImpl<int> to_int() const {
        return to_type<int>();
    }

    calcImpl<long long> to_long_long() const {
        return to_type<long long>();
    }

    calcImpl<unsigned int> to_unsigned_int() const {
        return to_type<unsigned int>();
    }

    calcImpl<unsigned long long> to_unsigned_long_long() const {
        return to_type<unsigned long long>();
    }

    calcImpl<double> to_double() const {
        return to_type<double>();
    }

    calcImpl<long double> to_long_double() const {
        return to_type<long double>();
    }

    calcImpl<size_t> to_size_t() const {
        return to_type<size_t>();
    }

    calcImpl<sf_count_t> to_sf_count_t() const {
        return to_type<sf_count_t>();
    }

    operator T() const {
        return value;
    }

private:
    bool is_integer() const {
        return std::is_integral<T>::value;
    }

    template <typename U>
    void check_overflow(T value) const {
        const long double ld_value = static_cast<long double>(value);
        if (ld_value > static_cast<long double>(std::numeric_limits<U>::max()) ||
            ld_value < static_cast<long double>(std::numeric_limits<U>::lowest())) {
            throw std::overflow_error("型への変換時にオーバーフローが発生しました");
        }
    }
};

template <typename T>
calcImpl<T> safeCalc(T v) {
    return calcImpl<T>(v);
}


template <typename T>
calcImpl<T> safeCalc(const calcImpl<T>& obj) {
    return calcImpl<T>(obj.value);
}

GPU費用の募金箱 ↓

ここから先は

0字

¥ 1,000

この記事が気に入ったらチップで応援してみませんか?