リアルタイムでマイク入力をフーリエ変換

Question

こちらのサイトを元にリアルタイムで録音するアプリケーションを作成したのですが、このマイク入力してる音声に対しFFTをかけてリアルタイムで周波数を取り出したいと考えています。

MM_WIM_DATAのメッセージを受け取ったときの処理にフーリエ変換を入れればよいのだろうとは思うのですが、bTempに一度録音データでいっぱいになったバッファをコピーしているので、このbTempにFFTをかけるのだと思っているのですが、BYTE型に対して数式の処理のFFTをかけるということがいまいちピンときていません。

それともFFTをかけるべき変数が違うのでしょうか。マイク入力を逐次FFTしていく方法を教えていただきたいです。

また、可能であれば、逐次的にFFTしたデータをCSVなどで出力する場合、録音を終了し、終了するまでに逐次的にかけていたFFTで得たデータの平均を出力すればよいのでしょうか。

#include <windows.h>
//#include <mmeapi.h>
#include "resource.h"

#pragma comment(lib, "Winmm.lib")

#define TITLE TEXT("REC")
#define SRATE 48000

LRESULT CALLBACK WndProc(HWND hWnd , UINT msg , WPARAM wp , LPARAM lp)
{
    //各種クラスの宣言
    static WAVEFORMATEX wfe;
    static HWAVEOUT hWaveOut;
    static HWAVEIN hWaveIn;
    static BYTE *bWave1 , *bWave2 , *bSave , *bTmp;
    static WAVEHDR whdr1 , whdr2;
	static DWORD dwLength = 0, dwCount;
	static BOOL blReset = FALSE;


    //メッセージ毎の挙動
    switch (msg)
    {
        //解放
        case WM_DESTROY:
            if(bSave) free(bSave);
            PostQuitMessage(0);
            return 0;
        //作成
        case WM_CREATE:
            //WAVEのフォーマットを定義
            //モノラル,48kHz,16bit
            wfe.wFormatTag = WAVE_FORMAT_PCM;
            wfe.nChannels = 1;
            wfe.nSamplesPerSec = SRATE;
            wfe.nAvgBytesPerSec = SRATE*2;  //48k * 1ch * 16bit /8
            wfe.wBitsPerSample = 16;
			wfe.nBlockAlign = 2;			//wfe.nChannels * wfe.wBitsPerSample / 8
            //書き込み、録音、再生はどれもnBlockAlignの整数倍で行う
			wfe.cbSize = 0;
            return 0;
        case WM_COMMAND:
            switch(LOWORD(wp))
            {
                case IDB_REC:
					bWave1 = (BYTE*)malloc(SRATE * 10);
					bWave2 = (BYTE*)malloc(SRATE * 10);

					whdr1.lpData = (char*)bWave1;
					whdr1.dwBufferLength = SRATE * 5;
					whdr1.dwBytesRecorded = 0;
					whdr1.dwFlags = 0;
					whdr1.dwLoops = 1;
					whdr1.lpNext = NULL;
					whdr1.dwUser = 0;
					whdr1.reserved = 0;

					whdr2.lpData = (char*)bWave2;
					whdr2.dwBufferLength = SRATE * 5;
					whdr2.dwBytesRecorded = 0;
					whdr2.dwFlags = 0;
					whdr2.dwLoops = 1;
					whdr2.lpNext = NULL;
					whdr2.dwUser = 0;
					whdr2.reserved = 0;

					waveInOpen(&hWaveIn, WAVE_MAPPER, &wfe, (DWORD)hWnd, 0, CALLBACK_WINDOW);
					waveInPrepareHeader(hWaveIn, &whdr1, sizeof(WAVEHDR));
					waveInPrepareHeader(hWaveIn, &whdr2, sizeof(WAVEHDR));
                    break;
                case IDB_PLAY:
					whdr1.lpData = (char*)bSave;
					whdr1.dwBufferLength = dwLength;
					whdr1.dwBytesRecorded = 0;
					whdr1.dwFlags = WHDR_BEGINLOOP | WHDR_ENDLOOP;
					whdr1.dwLoops = 1;
					whdr1.dwUser = 0;
					whdr1.reserved = 0;

                    waveOutOpen(&hWaveOut , WAVE_MAPPER , &wfe , (DWORD)hWnd , 0 , CALLBACK_WINDOW);
		        	waveOutPrepareHeader(hWaveOut , &whdr1 , sizeof(WAVEHDR));
					waveOutWrite(hWaveOut, &whdr1, sizeof(WAVEHDR));
		        	break;
				case IDB_REC_END:
					blReset = TRUE;
					waveInReset(hWaveIn);
					break;
				case IDB_PLY_END:
					waveOutReset(hWaveOut);
					break;
            }
        return 0;
        case MM_WIM_OPEN:
			dwLength = 0;
			bSave = (BYTE*)realloc(bSave, 1);

            EnableWindow(GetDlgItem(hWnd , IDB_PLAY) , FALSE);
            EnableWindow(GetDlgItem(hWnd , IDB_REC) , FALSE);
			EnableWindow(GetDlgItem(hWnd , IDB_REC_END) , TRUE);
            waveInAddBuffer(hWaveIn , &whdr1 , sizeof(WAVEHDR));
			waveInAddBuffer(hWaveIn , &whdr2, sizeof(WAVEHDR));
            waveInStart(hWaveIn);
            return 0;
        case MM_WIM_DATA:
			bTmp = (BYTE*)realloc(bSave, dwLength + ((LPWAVEHDR)lp)->dwBytesRecorded);
			if (blReset || !bTmp)
			{
				waveInClose(hWaveIn);
				blReset = FALSE;
				return 0;
			}
			bSave = bTmp;
			for(dwCount = 0; dwCount < ((LPWAVEHDR)lp)->dwBytesRecorded; dwCount++)
				*(bSave + dwLength + dwCount) = *(((LPWAVEHDR)lp)->lpData + dwCount);
			dwLength += ((LPWAVEHDR)lp)->dwBytesRecorded;
			waveInAddBuffer(hWaveIn, (LPWAVEHDR)lp, sizeof(WAVEHDR));
            return 0;
        case MM_WIM_CLOSE:
			waveInUnprepareHeader(hWaveIn, &whdr1, sizeof(WAVEHDR));
			waveInUnprepareHeader(hWaveIn, &whdr2, sizeof(WAVEHDR));
			free(bWave1);
			free(bWave2);

			EnableWindow(GetDlgItem(hWnd, IDB_PLAY), TRUE);
			EnableWindow(GetDlgItem(hWnd, IDB_REC), TRUE);
			EnableWindow(GetDlgItem(hWnd, IDB_REC_END), FALSE);
	    	return 0;
	   case MM_WOM_OPEN:
		   EnableWindow(GetDlgItem(hWnd, IDB_PLAY), FALSE);
		   EnableWindow(GetDlgItem(hWnd, IDB_REC), FALSE);
		   EnableWindow(GetDlgItem(hWnd, IDB_PLY_END), TRUE);
	    	return 0;
	    case MM_WOM_DONE:
	    	waveOutClose(hWaveOut);
	    	return 0;
    	case MM_WOM_CLOSE:
			waveOutUnprepareHeader(hWaveOut, &whdr1, sizeof(WAVEHDR));
			EnableWindow(GetDlgItem(hWnd, IDB_PLAY), TRUE);
			EnableWindow(GetDlgItem(hWnd, IDB_REC), TRUE);
			EnableWindow(GetDlgItem(hWnd, IDB_PLY_END), FALSE);
	    	return 0;
    }
    return DefWindowProc(hWnd , msg , wp , lp);

}

int WINAPI WinMain(HINSTANCE hInstance , HINSTANCE hPrevInstance , 
                    PSTR lpCmdLine , int nCmdShow)
{
    HWND hwnd;
    MSG msg;
    WNDCLASS winc;

    winc.style          = CS_HREDRAW | CS_VREDRAW;
    winc.lpfnWndProc    = WndProc;
    winc.cbClsExtra     = 0;
    winc.cbWndExtra     = DLGWINDOWEXTRA;
    winc.hInstance      = hInstance;
    winc.hIcon          = LoadIcon(NULL , IDI_APPLICATION);
    winc.hCursor        = LoadCursor(NULL , IDC_ARROW);
    winc.hbrBackground  = (HBRUSH)CreateSolidBrush(GetSysColor(COLOR_MENU));
    winc.lpszMenuName	= NULL;
	winc.lpszClassName	= TEXT("KITTY");

    if (!RegisterClass(&winc)) return -1;
    hwnd = CreateDialog(hInstance , TEXT("KITTY") , 0 , NULL);
    
    if (hwnd == NULL) return -1;

    while(GetMessage(&msg , NULL , 0 , 0)) {
        TranslateMessage(&msg);
        DispatchMessage(&msg);
    }
    return msg.wParam;

}

Answer 1

リアルタイムで周波数を取り出したいと考えています。

スペアナのようなものを想定しているのでしょうか？

でしたら、データの一部の時間区間を窓関数で切り抜いてフーリエ変換を行い、窓関数の時間軸をずらしながら上記の処理を繰り返していく『短時間フーリエ変換』という方法がありますので、こちらを使用するとよいのではないでしょうか。

Answer 2

    case MM_WIM_DATA:
        bTmp = (BYTE*)realloc(bSave, dwLength + ((LPWAVEHDR)lp)->dwBytesRecorded);
        if (blReset || !bTmp)
        {
            waveInClose(hWaveIn);
            blReset = FALSE;
            return 0;
        }
        bSave = bTmp;
        for(dwCount = 0; dwCount < ((LPWAVEHDR)lp)->dwBytesRecorded; dwCount++)
            *(bSave + dwLength + dwCount) = *(((LPWAVEHDR)lp)->lpData + dwCount);
        dwLength += ((LPWAVEHDR)lp)->dwBytesRecorded;
        waveInAddBuffer(hWaveIn, (LPWAVEHDR)lp, sizeof(WAVEHDR));
        return 0;

bTmpですと、0秒から現在までの録音できた全てのデータになってしまいますね。
普通リアルタイムでFFTするというとある一定時間をサンプリングしたものを用いますので
((LPWAVEHDR)lp)->lpDataが指す領域に((LPWAVEHDR)lp)->dwBytesRecordedbyte分格納されているデータをFFTするとよいと思います。

        //WAVEのフォーマットを定義
        //モノラル,48kHz,16bit

なのでwavの仕様から、BYTEよりもint16_tを使ったほうがよいでしょう。

Comments

[YMD_kts]

FFTで円周率を扱うのでdouble* が良いのではないかと思ったのですが、SHORTの方が、キレイに波形データのバイトとマッチしてて良いのですかね。

[asm]

最終および中間の結果がdoubleになるのはいいのですが
入力であるwavデータ自体は符号付き16bit整数が配列状になっているものであることに注意が必要です。

Answer 3

BYTE型に対して数式の処理のFFTをかけるということがいまいちピンときていません

そのピンと来ていないのがよくわからないのですが、FFTは通常2のn乗個のデータを渡して処理するものです。何個渡すかによって周波数の解像度がかわります。当然たくさん渡したほうが解像度はよくなりますが、それだけのデータが貯まるまで時間はかかるし、計算量も多くなります。

Windowsのそのへんの仕組みは知りませんが、bTempというのは何サンプルぐらいのデータが得られるのでしょう? FFTに渡したいデータ個数よりbTempのデータ数の方が少ない場合があるなら、別途バッファを設けてFFTにわたすデータを貯めてやらなきゃいけないかも知れません。