信号処理分野のための GNU Octave/MATLAB/Scilab 入門

---

本ページで公開している資料の対象者とコンセプト

このページには 『“信号処理分野のための GNU Octave/MATLAB/Scilab 入門”, 神戸市立工業高等専門学校 総合情報センター広報 第30号 (2018) pp.10-17』 関連のファイルおよびサポート情報を掲載しています。

• 少しはプログラミング経験があるが GNU Octave/MATLAB/Scilab は（ほぼ）初めて使う，という人を対象にしています。
• 工学的な話に特化しています。数学的な話は扱っていません。高専の電子工学科５年生を強く意識しています。
• 簡単で実用的なサンプルを集めているつもりですので，眺めるだけではなく実際に手を動かしてやってみて下さい。

MATLABの実行画面例

Scilabの実行画面例

解説資料本体のダウンロード

まずは資料本体を以下からダウンロードしてください。電子ブックリーダーやタブレット端末をお持ちの方は epub または mobi が便利です（画面サイズに応じてレイアウトされます）。

• フルバージョン PDF（8ページ） (602 kB)
• フルバージョン EPUB（電子ブック形式：iBooks, Sony Readerなど） (44 kB)
• フルバージョン mobi（電子ブック形式：Kindle） (181 kB)
• 抜粋バージョン PDF（4ページ） (446 kB)

インストール

以下のいずれか一つをインストールしてください。（MATLABをお持ちの方はもちろんそれでOK）

GNU Octave

GNU Octave のサイト で GNU Octave 4.4.1 を入手（有用であると感じたらプロジェクトに寄附することを勧める）

MATLAB

iOS/Android版のMATLAB Mobileでも無料で全く同じ機能が使える。特にタブレット端末ではある程度快適に利用できる。（クラウドでの実行なので端末は要オンライン。アカウント登録だけは必要だが無料なので登録しておいて損は無い。）
[ iOS版（無料） / Android版（無料） / 学生版の購入 / ホーム版の購入 (Signal Processing Toolbox必須) ]

Scilab

Scilab のサイト で Scilab 5.5.2 を入手（6.0.0でも構わない）（有用であると感じたらプロジェクトに寄附することを勧める）

課題の解答例（Octave/MATLAB対応）

答えの確認方法

点線の枠の中に

こんな感じで

答えが書かれています。マウスを持っていくと表示できます。
うまく表示されない場合は 30秒表示 ボタンを押して下さい。全ての答えが30秒間表示されます。
いきなり答えを見るのではなく，自身で考えてから見るようにして下さい。

◆課題１：dB→比 に戻してみる

例： 6 dB は 何倍か？

10.^(6/20)
（この例では . はなくても結果は同じです）

30秒表示

◆課題３：1から100までの合計は？

sum([1:100])
（この例では [ ] はなくても結果は同じです）

30秒表示

◆課題４：7～18までの整数それぞれの自乗の合計は？

sum([7:18].^2)

30秒表示

◆課題５５：-π ≦ x ＜ π の範囲で y = x と y = sin(x) の２つのグラフを重ねて表示してみる。（２つの関数の特徴が 明確になるように x の増分を充分細かく取る。また，凡例も付ける。）

x=[-pi:0.01:pi];
y1=x;
y2=sin(x);
plot(x,y1)
hold on
plot(x,y2)
grid on
legend('y = x','y = sin(x)')
xlim([-pi,pi])

30秒表示

◆課題６： x > 0 の範囲で y=sqrt(x), x ，x²，log(x) の４つのグラフを並べて表示してみる。（特徴を比較しやすいようにx の範囲と増分を調整する。また，各グラフにタイトルも付ける。）

x=[0.01:0.01:1.5];
y1=sqrt(x);
y2=x;
y3=x.^2;
y4=log(x);
subplot(2,2,1)
plot(x,y1)
grid on
title('y = √x')
subplot(2,2,2)
plot(x,y2)
grid on
title('y = x')
subplot(2,2,3)
plot(x,y3)
grid on
title('y = x^2')
subplot(2,2,4)
plot(x,y4)
grid on
title('y = log(x)')

% 行数を減らして見通しをよくするために
x=[0.01:0.01:1.5];
y1=sqrt(x);
y2=x;
y3=x.^2;
y4=x.^3;
subplot(2,2,1); plot(x,y1); grid on; title('y = √x');
subplot(2,2,2); plot(x,y2); grid on; title('y = x');
subplot(2,2,3); plot(x,y3); grid on; title('y = x^2');
subplot(2,2,4); plot(x,y4); grid on; title('y = log(x)');
% のように書くことも多い。
% （セミコロン ; を付けると並べることができる）

30秒表示

◆課題７：♪ミの音などを作成し，sig（♪ラ）と比べてみる。（sig は残しておく）

fs=8000;
T=0.2;
t=[0:T*fs-1]/fs;
f=440;
sig=sin(2*pi*f*t);
% ここまでは資料と同じ

f_mi=660;
sig_mi=sin(2*pi*f_mi*t);

plot(t,sig,'o-b'); hold on;
plot(t,sig_mi,'o-r');
xlabel('時刻 [s]'); ylabel('振幅 [arb. unit]');
legend('440 Hz','660 Hz')

sound(sig,fs);
sound(sig_mi,fs);

30秒表示

◆課題８：横軸をミリ秒にして課題７のグラフを描く（横軸の範囲も見やすく調整する：3.2.参照）。

plot(t*1000,sig,'o-b'); hold on;
plot(t*1000,sig_mi,'o-r');
xlim([0 20]); grid on;
xlabel('時刻 [ms]'); ylabel('振幅 [arb. unit]');
legend('440 Hz','660 Hz')

30秒表示

◆課題９：同じサンプリングレート（fs=8000Hz）で，3900Hzの信号ではどうなるか，グラフを描いて確認する。

fs=8000;
T=0.2;
t=[0:T*fs-1]/fs;
% ここまでは資料と同じ

f2=3900;
sig2=sin(2*pi*f2*t);
plot(t,sig2,'o-')
xlim([0 0.02]); grid on;
xlabel('時刻 [s]')
ylabel('振幅 [arb. unit]')
title(sprintf('%d Hz の正弦波を %d Hz でサンプリングした波形',f2,fs))

30秒表示

◆課題１０：sigの最大値，最小値，平均値，実効値をそれぞれ求めてみる。

fs=8000;
T=0.2;
t=[0:T*fs-1]/fs;
f=440;
sig=sin(2*pi*f*t);
% ここまでは資料と同じ

max(sig)
min(sig)
mean(sig)
sqrt(mean(sig.^2))

30秒表示

◆課題１１：一様分布乱数から疑似正規分布風の乱数を作る。

fs=8000;
r3=rand(1,fs)+rand(1,fs)+rand(1,fs)+rand(1,fs)+rand(1,fs);
mean(r3)
histogram(r3,20)
soundsc(r3,fs)

30秒表示

◆課題１２：各ファイルに保存された波形を Audacity / SoundEngine 等の波形編集ソフトで確認する。

・Audacity： Windows版ダウンロード / 公式サイト [Windows/Mac/Linux]

・SoundEngine Free： Windows版ダウンロード / 公式サイト [Windows]

◆課題１３：任意のwavファイルを読み込んで聞いてみる。また，fsを変えて試してみる。例えば，半音上げるにはどうすればよいか？

wavファイルが必要な場合：[モノラル (1.2MB) / ステレオ (2.4MB)]

[sig3,fs3]=audioread('mono.wav');
sound(sig3,fs3)
sound(sig3,fs3*(2.^(1/12)))

30秒表示

◆課題１４：4.3.の「sound(sig*2,fs)」では何が起きていたか？ グラフを書いて確認する（できればbeforeとafterを重ねて描いて比べる）。

fs=8000;
T=0.2;
t=[0:T*fs-1]/fs;
f=440;
sig=sin(2*pi*f*t);
% ここまでは資料と同じ

sig_clip=max(-1,min(1,sig*2));
plot(t,sig,'o-b');
hold on
plot(t,sig_clip,'o-r');
xlim([0 0.01]); ylim([-1.1 1.1]); grid on;
legend('before（オリジナル）','after（クリップ）')
xlabel('時刻 [s]')
ylabel('振幅 [arb. unit]')

sound(sig,fs)
sound(sig_clip,fs)

30秒表示

◆課題１５：4.1.で作成したsigのFFT演算をおこない，正しい横軸（点数ではなくHz）で表示する。

fs=8000;
T=0.2;
t=[0:T*fs-1]/fs;
f=440;
sig=sin(2*pi*f*t);
sig_clip=min(1,max(-0.8,min(0.8,sig)));
% ここまでは資料と同じ

% FFTの計算と横軸(周波数軸)の生成
sig_clip_fft=fft(sig_clip);
df=1/T;
freq_axis=[0:T*fs-1]*df;

% グラフ
stem(freq_axis,abs(sig_clip_fft),'o-');
title('クリップされた正弦波の振幅スペクトル')
xlabel('周波数 [Hz]')
ylabel('振幅スペクトル [arb. unit]')

30秒表示

◆課題１６：440Hzの正弦波を　(A) fs=8kHz，T=0.01s　(B) fs=10kHz, T=0.008s　の各条件で作成し，１つのグラフに表示する。縦軸のdB表記（最大値で正規化）も試す。また，このとき，それぞれのピーク値の周波数はどうなるか？

% (A)の信号の生成
fsa=8000;
Ta=0.01;
ta=[0:Ta*fsa-1]/fsa;
fa=440;
siga=sin(2*pi*fa*ta);

% (B)の信号の生成
fsb=10000;
Tb=0.008;
tb=[0:Tb*fsb-1]/fsb;
fb=440;
sigb=sin(2*pi*fb*tb);

% FFTの計算と横軸(周波数軸)の生成
siga_fft=fft(siga);
sigb_fft=fft(sigb);
dfa=1/Ta; freq_axis_a=[0:Ta*fsa-1]*dfa;
dfb=1/Tb; freq_axis_b=[0:Tb*fsb-1]*dfb;

% 波形の表示
subplot(3,1,1)
plot(ta,siga,'o-r'); hold on;
plot(tb,sigb,'o-b');
grid on; legend('(A)','(B)'); title('サンプリングされた波形');
xlabel('時刻 [s]'); ylabel('振幅 [arb. unit]')

% 振幅スペクトルのグラフ表示（リニア）
subplot(3,1,2);
plot(freq_axis_a,abs(siga_fft),'o-r'); hold on;   % 本来は stem で描くべきだが
plot(freq_axis_b,abs(sigb_fft),'o-b');               % 見にくいので plot で代用した
grid on; legend('(A)','(B)'); title('振幅スペクトル');
xlabel('周波数 [Hz]'); ylabel('振幅スペクトル [arb. unit]')

% 振幅スペクトルのグラフ表示（デシベル）（sigaのスペクトルの最大値で正規化）
siga_fft_normalized=siga_fft/max(abs(siga_fft));
sigb_fft_normalized=sigb_fft/max(abs(siga_fft));
subplot(3,1,3)
plot(freq_axis_a,20*log10(abs(siga_fft_normalized)),'o-r'); hold on;
plot(freq_axis_b,20*log10(abs(sigb_fft_normalized)),'o-b');
grid on; legend('(A)','(B)'); title('振幅スペクトル');
xlabel('周波数 [Hz]'); ylabel('正規化された振幅スペクトル [dB]');

30秒表示

◆課題１７：一様分布のノイズあるいは正規分布のノイズ（4.4.および4.5.参照）のスペクトルを調べてみる。また，インパルスのスペクトルも調べてみる。

fs=8000;
T=1;
t=[0:T*fs-1]/fs;
f=[0:T*fs-1]/T;

% 信号の生成
sig_noise_uniform = rand(1, T*fs) - 0.5;
sig_noise_normal = randn(1, T*fs);
sig_impulse = [1 zeros(1, T*fs-1)];

% スペクトルの計算
fft_noise_uniform = fft(sig_noise_uniform);
fft_noise_normal = fft(sig_noise_normal);
fft_impulse = fft(sig_impulse);

% 正規化スペクトル（dB）の計算
fft_noise_uniform_normalized = 20*log10(abs(fft_noise_uniform)/max(abs(fft_noise_uniform)));
fft_noise_normal_normalized = 20*log10(abs(fft_noise_normal)/max(abs(fft_noise_normal)));
fft_impulse_normalized = 20*log10(abs(fft_impulse)/max(abs(fft_impulse)));

% 波形表示
subplot(2,3,1); plot(t,sig_noise_uniform,'r');
ylim([-1 1]); title('波形（一様分布ノイズ）'); ylabel('振幅 [arb. unit]');
subplot(2,3,2); plot(t,sig_noise_normal,'b');
ylim([-5 5]); title('波形（正規分布ノイズ）'); xlabel('時刻 [s]');
subplot(2,3,3); plot(t,sig_impulse,'ko');
ylim([-1 2]); title('波形（インパルス）');

% 振幅スペクトル表示
subplot(2,3,4); plot(f,fft_noise_uniform_normalized,'r');
xlim([0 fs/2]); ylim([-60 6]); title('振幅スペクトル（一様分布ノイズ）'); ylabel('振幅スペクトル [dB]');
subplot(2,3,5); plot(f,fft_noise_normal_normalized,'b');
xlim([0 fs/2]); ylim([-60 6]); title('振幅スペクトル（正規分布ノイズ）'); xlabel('周波数 [Hz]')
subplot(2,3,6); plot(f,fft_impulse_normalized,'k');
xlim([0 fs/2]); ylim([-60 6]); title('振幅スペクトル（インパルス）');

30秒表示

◆課題１８：任意の音声や音楽に残響をたたみ込んで聞いてみる。

wavファイルが必要な場合：[モノラル音声 (1.2MB) / ホールの残響（インパルス応答）（ステレオ） (516kB)]

[sig_mono,fs_mono]=audioread('mono.wav');
[sig_ir,fs_ir]=audioread('ir_hall.wav');
if fs_mono~=fs_ir
    disp('サンプリングレートが異なりますので注意して下さい')
end

% 残響の左チャンネルだけを使ってモノラルの残響音を生成
sound(sig_mono,fs_mono)
sound(sig_ir,fs_ir)
sig_conv = conv(sig_mono,sig_ir(:,1));
sound(sig_conv,fs_mono)

% 以下，ステレオの残響音を生成する場合
sig_conv_L = conv(sig_mono,sig_ir(:,1));
sig_conv_R = conv(sig_mono,sig_ir(:,2));
sound([sig_conv_L sig_conv_R],fs_mono)

30秒表示

◆課題１９：推定した楽譜情報にノイズが含まれる理由は何か，考える。【ヒント：作成したスネアの音の特徴は？】

fs=8000;

% 1秒分のドラム演奏楽譜
gakufu = zeros(1,fs);
gakufu(round(0.2*fs))=1.0;
gakufu(round(0.4*fs))=0.5;
gakufu(round(0.5*fs))=0.5;
gakufu(round(0.6*fs))=2.0;
subplot(4,1,1); plot(gakufu); title('楽譜');

% スネアっぽい音の生成（ホワイトノイズが指数関数で減衰する）
SN=exp(-[0:fs-1]/500).*(rand(1,fs)-0.5);
subplot(4,1,2); plot(SN); title('スネア（単音）');
sound(SN, fs);

% 畳み込みによって出音を生成
deoto=conv(gakufu, SN);
subplot(4,1,3); plot(deoto); title('出音（演奏音）');
sound(deoto, fs);

% 相互相関によって演奏情報を推定
gakufu2=xcorr(deoto, SN);
subplot(4,1,4);plot(gakufu2);title('推定結果');

% ここまでは資料と同じ

% 元のスネア音の自己相関
SN_autocorr=xcorr(SN, SN);
figure
plot(SN_autocorr)
title('元のスネア音の自己相関（理想的なインパルスにはなっていない）')

30秒表示

◆課題２０：周波数 fc [Hz] の搬送波（正弦波）を周波数 fs [Hz] の信号で変調度 m で振幅変調する関数を作成し，波形と振幅スペクトルを確認する。

【 am.m の中身】
function [result,t] = am(fcar,fsig,m,fs,T);
% fcar: 搬送波（carrier）の周波数 [Hz]
% fsig: 信号波の周波数 [Hz]
% m : 変調度 [比]
% fs: サンプリングレート [Hz]
% T: 信号の時間長

t = [0:T*fs-1]/fs;
result = ( 1 + m * sin(2*pi*fsig*t) ) .* sin(2*pi*fcar*t);

【実行するプログラム】
fcar=440;
fsig=50;
m=0.7;
fs=8000;
T=0.5;

% 時間信号の生成と表示
[am_sig,t] = am(fcar,fsig,m,fs,T);
sound(am_sig,fs)
subplot(2,1,1);
plot(t,am_sig,'o-');
xlabel('時刻 [s]'); ylabel('振幅 [arb. unit]'); title('時間波形');

% 振幅スペクトルの計算と表示
f=[0:T*fs-1]/T;
am_sig_fft = fft(am_sig);
am_sig_fft_dB = 20*log10(abs(am_sig_fft) / max(abs(am_sig_fft)));
subplot(2,1,2);
plot(f,am_sig_fft_dB,'o-');
xlim([0 fs/2]); grid on;
xlabel('周波数 [Hz]'); ylabel('正規化された振幅 [dB]'); title('振幅スペクトル');

30秒表示

課題２１ 　 課題２２ 　 課題２３＆２４ 　 課題２５ 　 課題２６

さらに深く学びたい方へ：信号処理のサンプルプログラム

Octave/MATLAB/Scilab で動作する信号処理のデモプログラムを下記で公開している。いずれも実用的な信号処理手法を視覚的に理解できるように紹介しているので参考にして欲しい。なお，Octave/MATLAB 版は日本語版も用意している。
https://github.com/nagataniyoshiki/SignalProcessingDemo

著者について

長谷 芳樹 (博士（工学）） / Yoshiki NAGATANI (Ph.D.(Engineering))
神戸市立工業高等専門学校 電子工学科 准教授

  
   github.com/nagataniyoshiki/
   @nagataniyoshiki
   www.researchgate.net/profile/Yoshiki_Nagatani
   www.linkedin.com/in/nagataniyoshiki

  [ ウェブサイト: https://ultrasonics.jp/nagatani/ ]

---