d4cの出力にnanが含まれることがある

[y-chan]

こんにちは、音声合成の研究を行っているy-chanと申します。素晴らしいライブラリの提供をありがとうございます。

WorldのPythonラッパーであるPyWorldを用いて、音声分析を行なっていた際、d4cの出力結果である非周期性指標にnanが含まれる場合があることを確認しました。 これは、Worldでの再合成でも、別の手法を用いた場合でも再合成に大きな影響をもたらします。

しかし、nanを引き起こすwavファイルを、Worldのexample(examples/analysis_synthesis/analysis.cppをコンパイルしたバイナリ)で直接読み込むことはできませんでした。そのためexampleを改変し、double型の配列をバイナリから読み込むものをexamples/analysis_synthesis/analysis_dat.cppとして作成して調べた結果、World単体でも再現することを確認しました。

PyWorldでnanを引き起こしたファイルEMOTION100_010.wav.datを以下の手順で分析させたところ、出力結果にnanが含まれることがわかりました。

examples/analysis_synthesis/analysis_dat.cppを配置して、以下のコマンドを実行します。 当方の環境はMacですが、WindowsやLinuxでも同様に実行できるはずです。

cd World
mkdir build & cd build
cmake ..
make
cd example_bin
./analysis_synthesis_analysis_dat EMOTION100_010.wav.dat EM010.f0 EM010.sp EM010.ap

analysis_dat.cpp

このファイルでは、datファイルを読み込む以外にも、f0推定手法をHarvestに変更したり、`frame_period`を変更したりして、該当ファイルでnanが出る条件を合わせています。 また、doubleの配列バイナリとしてではなく、可読性の観点から数値を文字列として出力しています。 ```cpp //----------------------------------------------------------------------------- // Copyright 2016 seblemaguer // Author: https://github.com/seblemaguer // Last update: 2017/02/01 // // Summary: // The example analyzes a .wav file and outputs three files in each parameter. // Files are used to synthesize speech by "synthesis.cpp". // // How to use: // The format is shown in the line 253. //----------------------------------------------------------------------------- #include #include #include #include #include #include #include #if (defined (__WIN32__) || defined (_WIN32)) && !defined (__MINGW32__) #include #include #pragma comment(lib, "winmm.lib") #pragma warning(disable : 4996) #endif #if (defined (__linux__) || defined(__CYGWIN__) || defined(__MINGW32__) || defined(__APPLE__)) #include #include #endif // For .wav input/output functions. #include "audioio.h" // WORLD core functions. // Note: win.sln uses an option in Additional Include Directories. // To compile the program, the option "-I $(SolutionDir)..\src" was set. #include "world/d4c.h" #include "world/harvest.h" #include "world/matlabfunctions.h" #include "world/cheaptrick.h" #include "world/stonemask.h" #if (defined (__linux__) || defined(__CYGWIN__) || defined(__MINGW32__) || defined(__APPLE__)) // Linux porting section: implement timeGetTime() by gettimeofday(), #ifndef DWORD #define DWORD uint32_t #endif DWORD timeGetTime() { struct timeval tv; gettimeofday(&tv, NULL); DWORD ret = static_cast(tv.tv_usec / 1000 + tv.tv_sec * 1000); return ret; } #endif //----------------------------------------------------------------------------- // struct for WORLD // This struct is an option. // Users are NOT forced to use this struct. //----------------------------------------------------------------------------- typedef struct { double frame_period; int fs; double *f0; double *time_axis; int f0_length; double **spectrogram; double **aperiodicity; int fft_size; } WorldParameters; namespace { /** * Displaying wav file information * * @param fs : the sampling frequence * @param nbit : the number of bits to code the signal * @param x_length : the number of samples */ void DisplayInformation(int fs, int nbit, int x_length) { std::cout << "File information" << std::endl; std::cout << "Sampling : " << fs << " Hz " << nbit << " Bit" << std::endl; std::cout << "Length " << x_length << " [sample]" << std::endl; std::cout << "Length " << static_cast(x_length) / fs << "[sec]" << std::endl; } /** * F0 estimation function * * @param x : the signal samples * @param x_length : the number of samples * @param world_parameters : the world structure which is going to contains the F0 values (double format) * in the f0 structure field */ void F0Estimation(double *x, int x_length, WorldParameters *world_parameters) { HarvestOption option = {0}; InitializeHarvestOption(&option); // Modification of the option // When you You must set the same value. // If a different value is used, you may suffer a fatal error because of a // illegal memory access. option.frame_period = world_parameters->frame_period; // Valuable option.speed represents the ratio for downsampling. // The signal is downsampled to fs / speed Hz. // If you want to obtain the accurate result, speed should be set to 1. // option.speed = 1; // You should not set option.f0_floor to under world::kFloorF0. // If you want to analyze such low F0 speech, please change world::kFloorF0. // Processing speed may sacrify, provided that the FFT length changes. option.f0_floor = 71.0; option.f0_ceil = 1000.0; // You can give a positive real number as the threshold. // Most strict value is 0, but almost all results are counted as unvoiced. // The value from 0.02 to 0.2 would be reasonable. // option.allowed_range = 0.1; // Parameters setting and memory allocation. world_parameters->f0_length = GetSamplesForHarvest(world_parameters->fs, x_length, world_parameters->frame_period); world_parameters->f0 = new double[world_parameters->f0_length]; world_parameters->time_axis = new double[world_parameters->f0_length]; // double *refined_f0 = new double[world_parameters->f0_length]; std::cout << std::endl << "Analysis" << std::endl; DWORD elapsed_time = timeGetTime(); Harvest(x, x_length, world_parameters->fs, &option, world_parameters->time_axis, world_parameters->f0); std::cout << "Harvest: " << timeGetTime() - elapsed_time << " [msec]" << std::endl; // StoneMask is carried out to improve the estimation performance. // elapsed_time = timeGetTime(); // StoneMask(x, x_length, world_parameters->fs, world_parameters->time_axis, // world_parameters->f0, world_parameters->f0_length, refined_f0); // std::cout << "StoneMask: " << timeGetTime() - elapsed_time << " [msec]" << std::endl; // for (int i = 0; i < world_parameters->f0_length; ++i) // world_parameters->f0[i] = refined_f0[i]; // delete[] refined_f0; return; } /** * Spectral envelope estimation function * * @param x : the signal samples * @param x_length : the number of samples * @param world_parameters : the world structure which is going to contains the spectrogram (double values) * in the spectrogram structure field */ void SpectralEnvelopeEstimation(double *x, int x_length, WorldParameters *world_parameters) { CheapTrickOption option; InitializeCheapTrickOption(world_parameters->fs, &option); // This value may be better one for HMM speech synthesis. // Default value is -0.09. option.q1 = -0.15; // Important notice (2016/02/02) // You can control a parameter used for the lowest F0 in speech. // You must not set the f0_floor to 0. // It will cause a fatal error because fft_size indicates the infinity. // You must not change the f0_floor after memory allocation. // You should check the fft_size before excucing the analysis/synthesis. // The default value (71.0) is strongly recommended. // On the other hand, setting the lowest F0 of speech is a good choice // to reduce the fft_size. option.f0_floor = 71.0; // Parameters setting and memory allocation. world_parameters->fft_size = GetFFTSizeForCheapTrick(world_parameters->fs, &option); world_parameters->spectrogram = new double *[world_parameters->f0_length]; for (int i = 0; i < world_parameters->f0_length; ++i) { world_parameters->spectrogram[i] = new double[world_parameters->fft_size / 2 + 1]; } DWORD elapsed_time = timeGetTime(); CheapTrick(x, x_length, world_parameters->fs, world_parameters->time_axis, world_parameters->f0, world_parameters->f0_length, &option, world_parameters->spectrogram); std::cout << "CheapTrick: " << timeGetTime() - elapsed_time << " [msec]" << std::endl; } /** * Aperiodicity envelope estimation function * * @param x : the signal samples * @param x_length : the number of samples * @param world_parameters : the world structure which is going to contains the aperidicity (double values) * in the aperiodicity structure field */ void AperiodicityEstimation(double *x, int x_length, WorldParameters *world_parameters) { D4COption option; InitializeD4COption(&option); option.threshold = 0; // Parameters setting and memory allocation. world_parameters->aperiodicity = new double *[world_parameters->f0_length]; for (int i = 0; i < world_parameters->f0_length; ++i) { world_parameters->aperiodicity[i] = new double[world_parameters->fft_size / 2 + 1]; } DWORD elapsed_time = timeGetTime(); // option is not implemented in this version. This is for future update. // We can use "NULL" as the argument. D4C(x, x_length, world_parameters->fs, world_parameters->time_axis, world_parameters->f0, world_parameters->f0_length, world_parameters->fft_size, &option, world_parameters->aperiodicity); std::cout << "D4C: " << timeGetTime() - elapsed_time << " [msec]" << std::endl; } void DestroyMemory(WorldParameters *world_parameters) { delete[] world_parameters->time_axis; delete[] world_parameters->f0; for (int i = 0; i < world_parameters->f0_length; ++i) { delete[] world_parameters->spectrogram[i]; delete[] world_parameters->aperiodicity[i]; } delete[] world_parameters->spectrogram; delete[] world_parameters->aperiodicity; } } // namespace /** * Main function * */ int main(int argc, char *argv[]) { if (argc != 5) { std::cerr << argv[0] << " " << std::endl; return EXIT_FAILURE; } // 2016/01/28: Important modification. // Memory allocation is carried out in advanse. // This is for compatibility with C language. // int x_length = GetAudioLength(argv[1]); int x_length = 0; int fs, nbit; std::ifstream wav_dat(argv[1], std::ios::in | std::ios::binary ); if (!wav_dat) { std::cerr << "error: File \"" << argv[1] << "\" not found" << std::endl; return EXIT_FAILURE; } std::vector wav_list; while( true ) { double d = 0; // ファイルを読み込む wav_dat.read( (char *)&d, sizeof(double) ); if (wav_dat.eof()) { break; } // 読み込んだサイズだけ文字列へ追加 wav_list.push_back(d); x_length++; } // if (x_length < 0) { // std::cerr << "error: File \"" << argv[1] << "\" is not a .wav format" << std::endl; // return EXIT_FAILURE; // } double *x = wav_list.data(); fs = 24000; nbit = 16; // wavread() must be called after GetAudioLength(). // int fs, nbit; // wavread(argv[1], &fs, &nbit, x); DisplayInformation(fs, nbit, x_length); // 2016/02/02 // A new struct is introduced to implement safe program. WorldParameters world_parameters; // You must set fs and frame_period before analysis/synthesis. world_parameters.fs = fs; // 5.0 ms is the default value. // Generally, the inverse of the lowest F0 of speech is the best. // However, the more elapsed time is required. world_parameters.frame_period = 256 / (double)fs * 1000; //--------------------------------------------------------------------------- // Analysis part //--------------------------------------------------------------------------- // F0 estimation F0Estimation(x, x_length, &world_parameters); // Spectral envelope estimation SpectralEnvelopeEstimation(x, x_length, &world_parameters); // Aperiodicity estimation by D4C AperiodicityEstimation(x, x_length, &world_parameters); std::cout << "fft size = " << world_parameters.fft_size << std::endl; //--------------------------------------------------------------------------- // Saving part //--------------------------------------------------------------------------- // F0 saving std::ofstream out_f0(argv[2], std::ios::out | std::ios::binary); if(!out_f0) { std::cerr << "Cannot open file: " << argv[2] << std::endl; return EXIT_FAILURE; } out_f0.write(reinterpret_cast(world_parameters.f0), std::streamsize(world_parameters.f0_length * sizeof(double))); out_f0.close(); // Spectrogram saving std::ofstream out_spectrogram(argv[3], std::ios::out | std::ios::binary); if(!out_spectrogram) { std::cerr << "Cannot open file: " << argv[3] << std::endl; return EXIT_FAILURE; } // write the sampling frequency out_spectrogram.write(reinterpret_cast(&world_parameters.fs), std::streamsize( sizeof(world_parameters.fs) ) ); // write the frame period out_spectrogram.write(reinterpret_cast(&world_parameters.frame_period), std::streamsize( sizeof(world_parameters.frame_period) ) ); // write the spectrogram data for (int i=0; i(world_parameters.spectrogram[i]), std::streamsize((world_parameters.fft_size / 2 + 1) * sizeof(double))); } out_spectrogram.close(); // Aperiodicity saving // std::ofstream out_aperiodicity(argv[4], std::ios::out | std::ios::binary); std::ofstream out_aperiodicity(argv[4], std::ios::out); if(!out_aperiodicity) { std::cerr << "Cannot open file: " << argv[4] << std::endl; return EXIT_FAILURE; } for (int i=0; i

ここで、伺いたい点が2点あります。

1. この問題を根本的に修正する方法はありますでしょうか？
   • いくらか掘り下げてみたところ、中間表現のstatic_group_delayが怪しいように感じましたが、私には適切な修正方法がわかりませんでした。
2. 修正が難しい場合、d4cのnanの出力を任意の値で置き換えて使うしかないかと考えています。しかし、置き換えに適した値はわかりませんでした。置き換えに適した値は何かあるのでしょうか？
   • d4cの出力が、coarse aperiodicityを補完して生成しているものであるので、グラデーションのような値でうめるべきかと考えています。しかし、d4cの出力に対してアプローチすると実装が複雑化しそうで、coarse aperiodicityの時点でnanを消せると一番綺麗になるのかなと思っています。とはいえ、現時点でWorldはcoarse aperiodicityを出力していないので、アプローチが難しいのが現状です。

以上、ご回答いただければ幸いです。 よろしくお願いします。

[mmorise]

まだ解決はしていませんが，エラーは再現できてこれがバグか否かまでは確認できたので現状報告です． １．こちらはバグで修正可能と思われます． もともとの実装であるMATLAB版では再現しませんでしたので，C++へ移植した際に何らかの問題が生じているようです．

２，３への回答としてですが，暫定版としてcoarse_aperiodicityでNaNの部分を1.0に置き換えれば，MATLAB版に近い挙動になりそうです．

なんとなく原因は分かってきたのですが，まだ修正をどうするかまでが決まっておりません．修正法を含めて考えてみますので，少しお時間を頂きますようお願いいたします（おそらく何らかのセーフガードを入れることになります）．

[mmorise]

とりあえず，大まかな原因と解決法は分かったので共有いたします．簡単に言えば，入力信号にごく小さい（結果に影響を及ぼさないと思われる範囲でやや大きめな）雑音を加えればOKです．

for (int i = 0; i < x_length; ++i)
  x[i] = x[i] + 0.000001 * randn();

今回でしたら，上記のコードをAperiodicityEstimationの手前に入れれば動きました．

問題はそもそもの原因ですが，これは短時間フーリエ変換でパワースペクトルを求めた際に，ナイキスト周波数付近で極端に小さな値になっていることです．その帯域のパワーが途中の処理で完全な0になってしまっており，途中でパワーで割り算する処理があるためそこから破綻していくようです．そのため，上記のように全帯域にパワーを乗せることで回避そのものはできます．

元々この影響は想定しており，D4CでもGetWindowedWaveform内部関数で同様の処理を入れていましたが，ノイズ量が小さすぎて効果が発揮されなかったようです．これを増やすことでも解決できそうではありますが，この値は音声に雑音を加えることで根本的に非周期性成分を増やすことになりますので，エラーが生じないバランスをどうとるかについては議論の余地があります．

その辺についてはなんとも言いにくいのですが，もしご希望があればお知らせ頂ければ検討いたします．

[mmorise]

こちら，前のコメントで書いたとおり加算するノイズ量を調整することで解決することにしました．通常音声を対象に影響を調べたところ，概ね0.001%未満の誤差で収まりそうなので，推定結果に与える影響はないと思います．ノイズ量の定数は新たに追加したので，d4c.cppとconstantnumbers.hの2つを修正しました．こちらでは問題なくNaNが消えておりますので，ご確認いただけると助かります．

[y-chan]

返信遅くなりすみません！ 迅速な調査・修正助かります...! あまり知識が深くはないのですが、誤差が大きくならないのであれば、雑音成分を足すのが良いのかなと思いました。 後ほどこちらでも確認してみます。

ありがとうございます！！

[y-chan]

大変おそくなりました、こちらでも複数のファイルで試してみましたが、nanの出現がなくなったので、無事に解決できたものと思います。 こちらのIssueはcloseとしたいと思います！ ありがとうございました...!