LoginSignup
1
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@y__h

MATLABによる音声信号の周波数特性分析システムの実装と解説

Last updated at Posted at 2025-01-22

はじめに

音声信号の周波数分析は、音声処理や音響工学の分野で非常に重要な技術です。本記事では、MATLABを使用して音声信号の周波数特性を分析するシステムの実装方法と、その背後にある理論について詳しく解説します。

システムの概要

このシステムは、音声ファイル（.wav形式）を入力として受け取り、以下の処理を行います：

  1. 時間波形の表示と保存
  2. 周波数特性の分析と可視化
  3. 複数音声ファイルの一括処理
  4. 分析結果の自動保存

音声データの読み込みと基本構造

デジタル音声データの構造

音声データは、MATLABのaudioread関数を使用して読み込みます：

[y, Fs] = audioread(filePath);

この関数は2つの重要な要素を返します：

  1. y: 音声のサンプル値データ

    • 値の範囲は-1.0から1.0に正規化されています
    • モノラル音声の場合は1次元配列
    • ステレオ音声の場合は2列の行列（左右チャンネル）

  2. Fs: サンプリング周波数（Hz）

    • 一般的なCDクオリティの場合は44.1kHz
    • 1秒間に何回サンプリングするかを示す値

サンプリングと時間軸の生成

時間軸は以下のように生成します：

duration = info.Duration;  % 音声の長さ（秒）
t = linspace(0, duration, length(y));  % 時間軸の生成

周波数分析の実装

高速フーリエ変換（FFT）の適用

FFTは時間領域の信号を周波数領域に変換する数学的手法です。実装は以下のようになります：

N = length(y);  % 信号の長さ
Y = fft(y);     % FFTの実行

FFTの結果であるYは複素数の配列となり、以下の情報を含んでいます：

  • 各周波数成分の振幅情報（絶対値）
  • 位相情報（偏角）

周波数軸の生成

周波数軸は以下のように生成します：

f = Fs * (0:(N/2))/N;  % 周波数軸（0からFs/2まで）

ここでの重要なポイント：

  • 0Hzからナイキスト周波数（Fs/2）までの範囲をカバー
  • ナイキストの定理に基づく周波数の上限設定
  • 周波数分解能はFs/Nとなる

スペクトル計算と正規化

振幅スペクトルの計算は以下のように行います：

P2 = abs(Y/N);  % 振幅スペクトルの計算
P1 = P2(1:N/2+1);  % 片側スペクトルの取得
P1(2:end-1) = 2*P1(2:end-1);  % エネルギー補正

この処理で：

  • 複素数の絶対値を取得
  • 信号長で正規化
  • 片側スペクトルのエネルギーを保存

デシベル変換と表示

人間の聴覚特性に合わせて、対数スケール（dB）で表示します：

plot(f, 20*log10(P1), 'LineWidth', 1);
xlabel('Frequency (Hz)');
ylabel('Magnitude (dB)');
grid on;

システムの特徴と利点

マルチチャンネル対応

ステレオ音声の場合、左右チャンネルを個別に処理し、比較することができます：

if numChannels == 2
    Y1 = fft(y(:,1));  % 左チャンネル
    Y2 = fft(y(:,2));  % 右チャンネル
    % それぞれのチャンネルで同様の処理を実行
end

システムの実装

0. システムの全体像

全文表示 
% 音声分析・可視化プログラム
clear all; close all;

% 入出力ディレクトリの設定
inputDir = fullfile(pwd, 'database', 'input');
outputDir = fullfile(pwd, 'database', 'output');

% 入力ディレクトリから.wavファイルの一覧を取得
audioFiles = dir(fullfile(inputDir, '*.wav'));

% すべての波形と周波数特性を表示するための準備
numFiles = length(audioFiles);
figure('Position', [100 100 800 numFiles*200]); % ファイル数に応じて縦長の図を作成

% 各音声ファイルに対して処理を実行
for fileIdx = 1:length(audioFiles)
    % 音声ファイルの読み込み
    fileName = audioFiles(fileIdx).name;
    filePath = fullfile(inputDir, fileName);
    [y, Fs] = audioread(filePath);
    
    % 音声の基本情報を取得
    info = audioinfo(filePath);
    duration = info.Duration;
    numChannels = info.NumChannels;
    
    % このファイル用の出力ディレクトリを作成
    fileOutputDir = fullfile(outputDir, fileName(1:end-4));
    if ~exist(fileOutputDir, 'dir')
        mkdir(fileOutputDir);
    end
    
    % 時間軸の生成
    t = linspace(0, duration, length(y));
    
    % 時間波形の描画（個別ファイル用）
    figure('Position', [100 100 800 400]);
    if numChannels == 1
        plot(t, y, 'b', 'LineWidth', 1);
        ylabel('Amplitude');
    else
        % 各チャンネルの最大振幅を計算
        max_amp_1 = max(abs(y(:,1)));
        max_amp_2 = max(abs(y(:,2)));
        
        % 振幅の大きい順にプロット順序を決定
        if max_amp_1 < max_amp_2
            plot(t, y(:,2), 'Color', [1 0 0 0.4], 'LineWidth', 1.2);
            hold on;
            plot(t, y(:,1), 'Color', [0 0 1 0.4], 'LineWidth', 1.2);
            legend('Right Channel', 'Left Channel');
        else
            plot(t, y(:,1), 'Color', [0 0 1 0.4], 'LineWidth', 1.2);
            hold on;
            plot(t, y(:,2), 'Color', [1 0 0 0.4], 'LineWidth', 1.2);
            legend('Left Channel', 'Right Channel');
        end
        ylabel('Amplitude');
        grid minor;
        set(gca, 'GridAlpha', 0.2);
    end
    xlabel('Time (s)');
    grid on;
    
    % 波形プロットを保存
    waveformFile = fullfile(fileOutputDir, 'waveform.png');
    saveas(gcf, waveformFile);
    close gcf;
    
    % 周波数分析
    N = length(y);
    f = Fs * (0:(N/2))/N;
    
    % 周波数特性の描画（個別ファイル用）
    figure('Position', [100 100 800 400]);
    if numChannels == 1
        Y = fft(y);
        P2 = abs(Y/N);
        P1 = P2(1:N/2+1);
        P1(2:end-1) = 2*P1(2:end-1);
        plot(f, 20*log10(P1), 'b', 'LineWidth', 1);
    else
        Y1 = fft(y(:,1));
        Y2 = fft(y(:,2));
        P2_1 = abs(Y1/N);
        P2_2 = abs(Y2/N);
        P1_1 = P2_1(1:N/2+1);
        P1_2 = P2_2(1:N/2+1);
        P1_1(2:end-1) = 2*P1_1(2:end-1);
        P1_2(2:end-1) = 2*P1_2(2:end-1);
        plot(f, 20*log10(P1_1), 'b', f, 20*log10(P1_2), 'r');
        legend('Left Channel', 'Right Channel');
    end
    xlabel('Frequency (Hz)');
    ylabel('Magnitude (dB)');
    grid on;
    
    % 周波数特性プロットを保存
    spectrumFile = fullfile(fileOutputDir, 'spectrum.png');
    saveas(gcf, spectrumFile);
    close gcf;
    
    % 全ファイルの波形を1枚にまとめた図の作成
    figure(1);
    subplot(numFiles, 1, fileIdx);
    if numChannels == 1
        plot(t, y, 'b', 'LineWidth', 1);
    else
        if max_amp_1 < max_amp_2
            plot(t, y(:,2), 'Color', [1 0 0 0.4], 'LineWidth', 1.2);
            hold on;
            plot(t, y(:,1), 'Color', [0 0 1 0.4], 'LineWidth', 1.2);
            legend('Right Channel', 'Left Channel');
        else
            plot(t, y(:,1), 'Color', [0 0 1 0.4], 'LineWidth', 1.2);
            hold on;
            plot(t, y(:,2), 'Color', [1 0 0 0.4], 'LineWidth', 1.2);
            legend('Left Channel', 'Right Channel');
        end
    end
    grid on;
    grid minor;
    set(gca, 'GridAlpha', 0.2);
    if fileIdx == numFiles
        xlabel('Time (s)');
    end
    ylabel('Amplitude');
    
    % 処理状況を表示
    fprintf('Processed: %s\n', fileName);
end

% 波形まとめ図を保存
saveas(gcf, fullfile(outputDir, 'all_waveforms.png'));
close gcf;

% 全ファイルの周波数特性を1枚にまとめた図の作成
figure('Position', [100 100 800 numFiles*200]);
for fileIdx = 1:length(audioFiles)
    % 音声ファイルの再読み込み
    fileName = audioFiles(fileIdx).name;
    filePath = fullfile(inputDir, fileName);
    [y, Fs] = audioread(filePath);
    info = audioinfo(filePath);
    numChannels = info.NumChannels;
    
    % 周波数分析の再実行
    N = length(y);
    f = Fs * (0:(N/2))/N;
    
    subplot(numFiles, 1, fileIdx);
    if numChannels == 1
        Y = fft(y);
        P2 = abs(Y/N);
        P1 = P2(1:N/2+1);
        P1(2:end-1) = 2*P1(2:end-1);
        plot(f, 20*log10(P1), 'b', 'LineWidth', 1);
    else
        Y1 = fft(y(:,1));
        Y2 = fft(y(:,2));
        P2_1 = abs(Y1/N);
        P2_2 = abs(Y2/N);
        P1_1 = P2_1(1:N/2+1);
        P1_2 = P2_2(1:N/2+1);
        P1_1(2:end-1) = 2*P1_1(2:end-1);
        P1_2(2:end-1) = 2*P1_2(2:end-1);
        plot(f, 20*log10(P1_1), 'b', f, 20*log10(P1_2), 'r');
        legend('Left Channel', 'Right Channel');
    end
    grid on;
    if fileIdx == numFiles
        xlabel('Frequency (Hz)');
    end
    ylabel('Magnitude (dB)');
end

% 周波数特性まとめ図を保存
saveas(gcf, fullfile(outputDir, 'all_spectrums.png'));
close gcf;

fprintf('Analysis completed. Results saved in output directory.\n');

ディレクトリ構成配下のようにします:

.
├── database
│   ├── input
│   │   ├── sampleSound_1.wav
│   │   ├── sampleSound_2.wav
│   │   ├── sampleSound_3.wav
│   ├── output/
└── analyzeAudioSignal.m

1. システムの初期設定

システム全体の設定を行います：

% 音声分析・可視化プログラム
clear all; close all;

% 入出力ディレクトリの設定
inputDir = fullfile(pwd, 'database', 'input');
outputDir = fullfile(pwd, 'database', 'output');

% 入力ディレクトリから.wavファイルの一覧を取得
audioFiles = dir(fullfile(inputDir, '*.wav'));

このコードでは、作業環境をクリアにし、入出力ディレクトリを設定します。音声ファイルは'database/input'ディレクトリから読み込まれ、結果は'database/output'ディレクトリに保存されます。

2. 音声データの読み込みと時間軸の生成

各音声ファイルに対して以下の処理を行います：

% 音声ファイルの読み込みと基本情報の取得
fileName = audioFiles(fileIdx).name;
filePath = fullfile(inputDir, fileName);
[y, Fs] = audioread(filePath);
info = audioinfo(filePath);
duration = info.Duration;
numChannels = info.NumChannels;

% 時間軸の生成
t = linspace(0, duration, length(y));

時間軸の生成は、波形の可視化に重要な役割を果たします。linspace関数により、音声の開始時刻から終了時刻まで均等に分割された時間軸が生成されます。

3. 周波数分析の実装

周波数分析の核となるFFT（高速フーリエ変換）の処理を実装します：

% 周波数分析の準備
N = length(y);  % 信号の長さ
f = Fs * (0:(N/2))/N;  % 周波数軸の生成（0からFs/2まで）

% FFTの実行と周波数特性の計算
Y = fft(y);  % FFTの実行
P2 = abs(Y/N);  % 振幅スペクトルの計算
P1 = P2(1:N/2+1);  % 片側スペクトルの取得
P1(2:end-1) = 2*P1(2:end-1);  % エネルギー補正

このコードでは、以下の重要な処理が行われています：

  1. FFTによる時間領域から周波数領域への変換
  2. 振幅スペクトルの計算
  3. 片側スペクトルの取得とエネルギー補正

4. 可視化と結果の保存

分析結果を視覚的に表現し、保存します：

% 周波数特性の描画
figure('Position', [100 100 800 400]);
plot(f, 20*log10(P1), 'b', 'LineWidth', 1);
xlabel('Frequency (Hz)');
ylabel('Magnitude (dB)');
grid on;

% 結果の保存
spectrumFile = fullfile(fileOutputDir, 'spectrum.png');
saveas(gcf, spectrumFile);

周波数特性は対数スケール（dB）で表示され、人間の聴覚特性により近い形で可視化されます。

まとめ

本システムは、音声信号の周波数特性を詳細に分析し、視覚化するための総合的なソリューションを提供します。FFTを核とする周波数分析の理論と実装を理解することで、より高度な音声処理アプリケーションの開発にも応用できます。

システムの拡張性も高く、必要に応じて以下のような機能を追加することも可能です：

  • より詳細な統計分析
  • リアルタイム処理への対応
  • 異なる窓関数の適用
  • スペクトログラム表示

音声信号処理の基礎となる本システムを理解し、活用することで、さまざまな音響分析のニーズに対応することができます。

参考文献・リソース

1
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
1
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?