[ TOPへ ]

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


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

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


MATLABの実行画面例
MATLABの実行画面例
Scilabの実行画面例
Scilabの実行画面例

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

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

インストール

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

GNU Octave
GNU Octave のサイト で GNU Octave 4.2.1 を入手(有用であると感じたらプロジェクトに寄附することを勧める)
MATLAB
iOS/Android版のMATLAB Mobileでも無料で全く同じ機能が使える。特にタブレット端末ではある程度快適に利用できる。(クラウドでの実行なので端末は要オンライン。アカウント登録だけは必要だが無料なので登録しておいて損は無い。)
[ iOS版(無料) / Android版(無料) / 学生版の購入 / ホーム版の購入 (Signal Processing Toolbox必須) ]
Scilab
Scilab のサイト で Scilab 5.5.2 を入手(6.0.0でも構わない)(有用であると感じたらプロジェクトに寄附することを勧める)

課題の解答例(Octave/MATLAB対応)

答えの確認方法

点線の枠の中に

こんな感じで

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

◆課題1:dB→比 に戻してみる

例: 6 dB は 何倍か?

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

◆課題3:1から100までの合計は?

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

◆課題4:7~18までの整数それぞれの自乗の合計は?

sum([7:18].^2)

◆課題55:-π ≦ x < π の範囲で y = x と y = sin(x) の2つのグラフを重ねて表示してみる。(2つの関数の特徴が 明確になるように 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])

◆課題6: x > 0 の範囲で y=sqrt(x), x ,x2,log(x) の4つのグラフを並べて表示してみる。(特徴を比較しやすいように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)');
% のように書くことも多い。
% (セミコロン ; を付けると並べることができる)

◆課題7:♪ミの音などを作成し,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);

◆課題8:横軸をミリ秒にして課題7のグラフを描く(横軸の範囲も見やすく調整する: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')

◆課題9:同じサンプリングレート(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))

◆課題10: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))

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

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)

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

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

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

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

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

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

◆課題14: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)

◆課題15: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]')

◆課題16:440Hzの正弦波を (A) fs=8kHz,T=0.01s (B) fs=10kHz, T=0.008s の各条件で作成し,1つのグラフに表示する。縦軸の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]');

◆課題17:一様分布のノイズあるいは正規分布のノイズ(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('振幅スペクトル(インパルス)');

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

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)

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

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('元のスネア音の自己相関(理想的なインパルスにはなっていない)')

◆課題20:周波数 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('振幅スペクトル');

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

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

著者について

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

  email e-mail
  GitHub github.com/nagataniyoshiki/
  Twitter @nagataniyoshiki
  ResearchGate www.researchgate.net/profile/Yoshiki_Nagatani
   LinkedIn www.linkedin.com/in/nagataniyoshiki

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