データI/Oと解析編 - ファイルや外部機器へのアクセスとその解析(2)

前回最後に示したモデルを実行するとシリアルデータを受信し、Time Scopeブロックに入力した信号の波形が表示されます。モータ位置の目標値入力がSIN波だったので、モータ軸の回転角をポテンショメータの電圧値で取得しているArduino UnoのADコンバータ値もSIN波となっています。Time Scopeブロックには波形の解析機能がありますので、それを利用して周波数を求めてみましょう。モデルで設定した通り約5秒の周期で動いていたようです。

モータの挙動に何Hzの信号成分が含まれているのか、不要な高調波成分を持たないか周波数特性解析を行ってみましょう。

tmp = logsout.get(1)            % Simulinkデータの取得
len = numel(tmp.Values.Data)    % データサイズ取得
data = double(squeeze(tmp.Values.Data(len-2047:end)));  % データ型やサイズの変更
data = detrend(data, 'constant');   % DCオフセットの除去

[sdata1, f1] = periodogram(data, hann(2048), 2048, Fs);     % ピリオドグラムの計算

figure
plot(f1, 20*log10(sdata1),'b'),grid on      % 周波数特性の表示
title('Power Spectrum Density'),xlabel('Freq (Hz) '),ylabel('Power (dB)')

1Hz未満のところでピークが立っているのが見えます。この周波数が何Hzなのか求めます。求めた周波数はマーカーと共に元の波形に重ね書きします。

[Pks,Locs] = findpeaks(20*log10(sdata1),'NPeaks',5,'Sortstr','descend'); % ピーク検出
peak_freq = f1(Locs)    % ピーク周波数の表示
hold on
plot(peak_freq,Pks,'ro')    % ピーク周波数にマーカー表示
text(peak_freq,Pks,num2str(peak_freq))  % 周波数文字列を表示
hold off

0.195Hzつまり約5秒周期で、元モデルで設定した通りに動作していることが確認できました。

Arduino Unoに搭載されているA/Dコンバータのサンプリングレートはそれほど高速ではないため高い周波数の信号を解析するのは困難です。高速サンプリングができるデータ収集ボードを使ってMATLABにデータを取得してみましょう。オプション製品Data Acquisition Toolboxを使うとデータ収集ボードからアナログやディジタル信号の入出力を行うことができます。今回はDigilentが販売するAnalog Discoveryという比較的安価な入門者用のデータ収集ボードを利用しました。このデバイスは、アナログ入出力が最大100MSPS、帯域5MHzで2ch、ディジタル入力も最大100MSPSで16chのデータを取得できます。

Analog DiscoveryをUSBで接続し、MATLABでセットアップを行います。

>> daq.getDevices       % 接続されているデバイス情報を表示
ans = 
digilent: Digilent Inc. Analog Discovery Kit Rev. C (Device ID: 'AD1')
   Analog input subsystem supports:
      -2.5 to +2.5 Volts,-25 to +25 Volts ranges
      Rates from 0.1 to 1000000.0 scans/sec
      2 channels ('1','2')
      'Voltage' measurement type

   Analog output subsystem supports:
      -5.0 to +5.0 Volts range
      Rates from 0.0 to 1000000.0 scans/sec
      2 channels ('1','2')
      'Voltage' measurement type


Properties, Methods, Events

>> sDaq = daq.createSession('digilent') % セッションの定義
sDaq = 
Data acquisition session using Digilent Inc. hardware:
   Will run for 1 second (10000 scans) at 10000 scans/second.
   No channels have been added.

Properties, Methods, Events

>> Fs_dq = 100e3;      % サンプリング周波数100kHz
>> ch = sDaq.addAnalogInputChannel('ad1', 1, 'Voltage');   % 入力チャンネル追加
>> ch = sDaq.addAnalogOutputChannel('ad1', 1, 'Voltage'); % 出力チャンネル追加
>> sDaq.Rate = Fs_dq;  % サンプリング周波数の設定
>> sDaq.Channels(1).Range = [-2.5 2.5];    % 入力電圧レンジ設定
>> sDaq.Channels(2).Range = [-5 5];       % 出力電圧レンジの設定

動作確認のためにアナログ入出力を直接接続して、アナログ出力したスイープ信号を入力から取得します。

t= 0:1/Fs_dq:1;                     % 1kHz のサンプルレートで 1 秒間
fo=25;f1=50e3;                     % 25Hz から始まり、50kHz まで増加
outData = 2*chirp(t,fo,1,f1,'q',[],'logarithmic')';     % スイープ信号生成
sDaq.queueOutputData(outData);          % データのキュー 
[data, timestamps, triggerTime] = sDaq.startForeground;     % 取得開始

取得した信号を時間軸応答とスペクトログラムで表示します。

figure
plot(timestamps, data);
figure
plot(timestamps, data), grid on
xlabel('Time (seconds)'); ylabel('Voltage (Volts)');
title('Sweep Signal Acquired with Analog Discovery')
xlim([0.06 0.08]), ylim([-2 2])  % 表示範囲の設定

figure
spectrogram(data,512,200,512,Fs_dq)     % スペクトログラムの表示

スペクトログラムは短時間でFFTした信号を時系列に並べたもので、周波数の時間変化を解析するために使用される手法です。上図では、横軸が周波数、縦軸が時間、表示色は信号の振幅を示しています。信号の振幅はグラフ右にあるカラーバーを参照して読み取ることが出来ます。プログラムしたとおり、Sin波のスイープ信号が出力されているのが確認できました。

次は音声信号の取り込み例です。ほとんどのWindows PCにはサウンドカードが内蔵されていますので、そのマイク端子またはライン端子からモノラル/ステレオ信号の取得をしてみましょう。

使用可能な機器としては、PC内蔵サウンドカード以外にも、USBやIEEE1394接続で使用するような汎用オーディオインタフェースを使用することができます。次のようにコマンド入力すると音声データを取得することが出来ます。

>> recObj = audiorecorder(44100, 16, 2)   % 44.1kHz 16bit 2chのオブジェクト作成
>> recordblocking(recObj, 5);
>> play(recObj);   % 取得したデータを再生
>> myRecording = getaudiodata(recObj);     % データ取得
>> figure, plot(myRecording);  % Plot

MATLAB基本機能では、Windows標準ドライバWindows Direct Soundを通してモノラル/ステレオの2chまでの信号を取得することが出来ます。さらにオプション製品DSP System Toolboxを使うと、次のように高機能な音声取り込み機能を利用することが出来ます。

• 音楽制作で使用されるASIOやWDM-KSなどの高性能なドライバが使用可能
• 3ch以上のマルチチャンネルで信号の取得
• ストリーミング処理(入出力をしながら例えばフィルタなどの処理が可能)

沢山のマイクを並べたマイクロホンアレイで信号を取得しながら処理や解析を行うといったことが出来るようになります。

DSP System Toolboxの機能を使ったプログラム例を見てみましょう。

% オーディオレコーダのオブジェクト作成
hain = dsp.AudioRecorder('SampleRate', 44100, 'NumChannels',2, 'DeviceName', 'UA-25EX'); 
% FIRフィルタのオブジェクト作成
hf = dsp.DigitalFilter('TransferFunction','FIR (all zeros)', 'Numerator', fir1(20,0.2))
% TimeScopeオブジェクト作成
hts = dsp.TimeScope('SampleRate', 44100,'ShowGrid', true);

tic
while toc < 10  % 10秒間実行
    out = step(hain);       % Audio In
    out = step(hf, out);     % Filtering
    step(hts,out);         % Plot
end

コマンドdsp.AudioRecorderは、PCに接続されているオーディオデバイスから信号を取得するためのオブジェクトを定義しています。同様にdsp.DigitalFilterとdsp.TimeScopeコマンドはディジタルフィルタと表示用のTimeScopeのオブジェクトを定義しています。TimeScopeの表示画面は下図に示しています。次にwhileループの中にあるstepコマンドで、定義したオブジェクトを繰返し実行して、ストリーミング処理により一連の処理である信号取得⇒処理⇒表示を行っています。dsp.TimeScopeは表示した信号の統計的な解析を行う機能を持っており、ピーク値、RMS、立ち上がり時間などを測定することが出来ます。

DSP System Toolboxでは、Simulinkでも音声信号の入出力を行ったり、外部ハードウェアデバイスでパラメータ調整を行ったりするための機能が提供されています。Simulinkモデルの例を見てみましょう。このモデルは、MIDIコントローラの制御信号を入力として、Simulinkブロックを使用した演算処理によりシンセサイザーのように信号を発振させ、オーディオデバイスから音声出力を行うようになっています。

動作の様子をビデオでご覧下さい。MIDIコントローラは音楽制作などに使われる汎用のインターェースで、そのコントロール信号をSimulinkデータとして出力することが出来ます。また、Simulink上で音声出力を行う機能としては、To Audio Deviceブロックが提供されています。他にも音声入力用のFrom Audio Deviceブロック、UDP Send/Receiveブロックなどが提供されており、外部ハードウェアを利用したデータ入出力に役立ちます。

著者紹介

松本 充史(まつもと　あつし)
MathWorks Japan
アプリケーションエンジニアリング部
シニアアプリケーションエンジニア
Mathworks JapanではMATLABの中でも特に信号処理やコード生成に関する機能を担当している。

今回の記事で紹介しているMATLABコードやSimulinkモデルのファイルをご要望の方はこちらまでお問い合わせ下さい。