MATLAB: SOFA形式のHRTFを用いた音声信号への畳み込み処理の実装方法

# 1. はじめに

頭部伝達関数（Head-Related Transfer Function: HRTF）は、音源から両耳までの音の伝達特性を表す関数です。人間は、この伝達特性の違いを利用して音源の方向を知覚しています。

HRTFを音声信号に適用することで、2チャンネルのステレオ音声でも、あたかも特定の方向から音が聞こえてくるような効果を作り出すことができます。この技術は、バーチャルリアリティ（VR）や3Dオーディオなど、様々な分野で活用されています。

本記事では、MATLABを使用してHRTFを音声に適用する方法を解説します。具体的には、SOFA形式のHRTFデータを読み込み、任意の方向から音が聞こえるように音声を加工するプログラムの実装方法を紹介します。

**対象読者**
- MATLABの基本的な使い方を理解している方
- 音響信号処理に興味がある方
- HRTFを使った音声の空間化を実装したい方

# 2. 開発環境の準備

## 2.1 SOFAツールボックスのインストール

まず、HRTFデータを扱うために必要な **SOFA（Spatially Oriented Format for Acoustics）** ツールボックスをインストールします。

1. SOFAツールボックスのダウンロード
    - [SOFAの公式GitHubページ](https://github.com/sofacoustics/SOFAtoolbox)にアクセス
    - `Downloads` チャプターにある **Current releases of SOFA Toolbox can be found on SourceForge** の `SourceForge` リンクをクリック
    
    ![SOFA_github.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3638255/1dca5960-cdb6-b023-9731-54591e4b4838.png)

    - SourceForgeページの黄緑色の `Download` ボタンをクリック
    
    ![SourceForge.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3638255/830c875f-27f7-0bde-f259-52469c439643.png)

    - 数秒待つとダウンロードが自動的に開始
    - `SOFA Toolbox for Matlab and Octave (バージョン)` というzipファイルを解凍

2. インストール手順
    - 解凍したファイルを適切なディレクトリに配置
      - 推奨: `C:\Users\(userName)\MATLAB` などのユーザーディレクトリ ( `MATLAB`フォルダは任意で作成して下さい)
    :::note warn
    注意
    `C:\Program Files\MATLAB\...\toolbox` への配置は権限の問題でエラーが発生する可能性があります
        :::
    - MATLABを起動
    - **ホーム** タブの **環境** ビューから `パスの設定` を選択
    
    ![matlab.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3638255/980634aa-6f83-270f-2e6b-0f4927625aed.png)

    - `フォルダを追加...` から解凍したフォルダ内の **SOFAtoolbox** フォルダを選択
    ![パスの設定後.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3638255/9f56ccb2-b9b1-4a9c-f52c-8360ca119525.png)

    - `保存` をクリックしてパスを確定

## 2.2 HRTFデータの入手

本実装では、東北大学が公開しているHRTFデータセットを使用します。

1. データのダウンロード
    - [東北大学のHRTFデータベース](https://www.riec.tohoku.ac.jp/pub/hrtf/index.html)にアクセス
    - **HRTF Data..........Download file of HRTF data and anthropometric data** の `HRTF Data` をクリック
    
    ![The RIEC HRTF Dataset.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3638255/4055cd63-761b-bec2-4098-4576cd77c940.png)

    - 任意のsubjectのHRTFファイルを選択（例：sub001）
    
    ![Download HRTF.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3638255/81b916e7-d529-2448-286a-845174b6c807.png)


**東北大学以外のHRTFデータベースを使用する場合:**
より詳細なHRTFデータベースの情報は[立体音響研究のための公開頭部伝達関数（HRTF）データベースとその動向](https://www.jstage.jst.go.jp/article/jasj/78/8/78_449/_pdf)を参照してください。

# 3. プログラム構造の解説

## 3.1 ディレクトリ構成

プロジェクトは以下のような階層構造で構築されています：

```planetext
.
├── database
│   ├── hrtf
│   │   └── RIEC_hrir_subject_001.sofa
│   ├── input
│   │   └── sample_sound.wav
│   └── output
├── str
│   ├── Config.m
│   ├── AudioProcessor.m
│   ├── HrtfManager.m
│   └── FileManager.m
└── main.m
```

各ディレクトリの役割：
- **database**: データファイルの格納
  - **hrtf**: HRTFファイルの保存
  - **input**: 入力音声ファイルの保存
  - **output**: 処理結果の出力先
- **str**: プログラムのソースファイル
- **main.m**: メインプログラム

## 3.2 クラス構成

プログラムは以下の4つのクラスで構成されています：

1. **Config**
    - 設定値の一元管理
    - ディレクトリパスやHRTF設定の定数定義
    ```matlab
    classdef Config
        properties (Constant)
            % ディレクトリ設定
            INPUT_DIR = fullfile(pwd, 'database', 'input');
            OUTPUT_DIR = fullfile(pwd, 'database', 'output');
            HRTF_DIR = fullfile(pwd, 'database', 'hrtf');
    
            % HRTF設定
            HRTF_FILE = 'RIEC_hrir_subject_001.sofa';
    
            % 角度設定
            START_ANGLE = 0;
            END_ANGLE = 360;
            ANGLE_STEP = 5;
        end
    end
    ```

2. **AudioProcessor**
    - 音声処理の実装
    - モノラル変換やHRTF畳み込み処理
    ```matlab
    classdef AudioProcessor
        methods (Static)
            function outputAudio = processAudio(audio, hrtfL, hrtfR)
                % 入力検証
                validateattributes(audio, {'numeric'}, {'2d'});
                validateattributes(hrtfL, {'numeric'}, {'vector'});
                validateattributes(hrtfR, {'numeric'}, {'vector'});
    
                % HRTFを畳み込み
                outputLeft = conv(audio, hrtfL);
                outputRight = conv(audio, hrtfR);
    
                % ステレオに変換して出力
                outputAudio = [outputLeft, outputRight];
            end
    
            function monoAudio = convertToMono(audio)
                monoAudio = mean(audio, 2);
            end
        end
    end
    ```

3. **HrtfManager**
    - HRTFデータの読み込みと管理
    - 指定角度のHRTFデータの取得
    ```matlab
    classdef HrtfManager
        methods (Static)
            function hrtfData = loadHrtfData()
                sofaFilePath = fullfile(Config.HRTF_DIR, Config.HRTF_FILE);
                if ~exist(sofaFilePath, 'file')
                    error('HRTFファイルが見つかりません: %s', sofaFilePath);
                end
                hrtfData = SOFAload(sofaFilePath);
            end
    
            function [hrtfL, hrtfR] = getHrtfData(hrtfData, azimuth, elevation)
                idx = find(hrtfData.SourcePosition(:,1)==azimuth & ...
                          hrtfData.SourcePosition(:,2)==elevation);
    
                if isempty(idx)
                    error('HRTF方向が見つかりません（方位角：%d度、仰角：%d度）', ...
                          azimuth, elevation);
                end
    
                hrtfL = reshape(squeeze(hrtfData.Data.IR(idx,1,:)), [], 1);
                hrtfR = reshape(squeeze(hrtfData.Data.IR(idx,2,:)), [], 1);
    
                if ~isnumeric(hrtfL) || ~isnumeric(hrtfR)
                    error('HRTFデータが不正です');
                end
            end
        end
    end
    ```

4. **FileManager**
    - ファイル操作の実装
    - 入出力ファイルの管理
    ```matlab
    classdef FileManager
        methods (Static)
            function wavFiles = getWavFiles()
                wavFiles = dir(fullfile(Config.INPUT_DIR, '*.wav'));
                if isempty(wavFiles)
                    error('入力フォルダにWAVファイルが見つかりません');
                end
            end
    
            function createOutputDir(outputDirPath)
                if ~exist(outputDirPath, 'dir')
                    mkdir(outputDirPath);
                end
            end
    
            function saveAudio(filepath, audio, fs)
                audiowrite(filepath, audio, fs);
            end
        end
    end
    ```

## 3.3 メインプログラム

メインプログラムの処理フロー：

1. 環境の初期化
2. HRTFデータの読み込み
3. 入力音声ファイルの処理
4. 各方向でのHRTF適用
5. 結果の保存

```matlab
function main()
    % パスの設定
    addpath(fullfile(pwd, 'str'));

    try
        % SOFAライブラリの初期化
        SOFAstart;

        % HRTFデータの読み込み
        hrtfData = HrtfManager.loadHrtfData();

        % 音声ファイルの処理
        processWavFiles(hrtfData);

    catch ME
        fprintf('エラーが発生しました: %s\\n', ME.message);
        return;
    end
end

function processWavFiles(hrtfData)
    wavFiles = FileManager.getWavFiles();

    for i = 1:length(wavFiles)
        processWavFile(wavFiles(i), hrtfData);
    end
end

function processWavFile(wavFile, hrtfData)
    filename = wavFile.name;
    filePath = fullfile(Config.INPUT_DIR, filename);

    % 音声ファイルの読み込みと処理
    [audio, fs] = audioread(filePath);
    audio = AudioProcessor.convertToMono(audio);

    % 出力ディレクトリの作成
    outputDir = fullfile(Config.OUTPUT_DIR, filename(1:end-4));
    FileManager.createOutputDir(outputDir);

    % 各角度での処理
    for angle = Config.START_ANGLE:Config.ANGLE_STEP:(Config.END_ANGLE-Config.ANGLE_STEP)
        processAngle(audio, fs, angle, hrtfData, outputDir, filename);
    end

    fprintf('音声ファイル %s の処理が完了しました。\\n', filename);
end

function processAngle(audio, fs, angle, hrtfData, outputDir, filename)
    [hrtfL, hrtfR] = HrtfManager.getHrtfData(hrtfData, angle, 0);
    outputAudio = AudioProcessor.processAudio(audio, hrtfL, hrtfR);

    outputFilename = sprintf('%s_%03d.wav', filename(1:end-4), angle);
    outputPath = fullfile(outputDir, outputFilename);
    FileManager.saveAudio(outputPath, outputAudio, fs);
end
```

# 5. 使用方法とカスタマイズ

## 5.1 基本的な使用方法

1. 入力音声ファイルを `database/input` に配置
2. `Config.m` で必要なパラメータを設定
3. `main.m` を実行
4. `database/output` から結果を取得

## 5.2 カスタマイズのポイント

- **角度設定の変更**
  - `Config.m` の `START_ANGLE`、`END_ANGLE`、`ANGLE_STEP` を調整
- **異なるHRTFの使用**
  - `Config.m` の `HRTF_FILE` を変更
  - 必要に応じて `HrtfManager.m` の読み込み処理を調整

# 参考文献・リソース

1. [SOFA - Spatially Oriented Format for Acoustics](https://github.com/sofacoustics/SOFAtoolbox)

2. [東北大学 RIEC HRTF Dataset](https://www.riec.tohoku.ac.jp/pub/hrtf/index.html)

3. [立体音響研究のための公開頭部伝達関数（HRTF）データベースとその動向](https://www.jstage.jst.go.jp/article/jasj/78/8/78_449/_pdf)

4. [MATLABでHRTF～頭部伝達関数とは？～聴覚の仕組みは解明されていない～](https://note.com/leftbank/n/n83db3a4d6108#828d24e5-5de9-400c-a1cc-ee2b4a1e1e73)