Contents

% 　フーリエ変換の練習,   aito, 2007/10/3

clear all; clc;

元信号の作成

t = 0 : 1/44100 : 5;
y = 0.4 * sin(2 * pi * 550 * t) + 0.3 * sin(2 * pi * 1100 * t) + 0.2 * sin(2 * pi * 1650 * t);
figure(1);
plot(t,y);
xlim([0, 0.01]);
title('SOURCE DATA');
xlabel('TIME [s]');
ylabel('INTENSITY [a.u.]');
wavwrite(y, 44100, 16, 'test2.wav');

フーリエ変換と窓関数を用いた短時間フーリエ変換との比較

[y, fs, bits] = wavread('test2.wav'); % sample1.wav を読み込む
Y = fft(y)/length(y); % フーリエ変換
Y = Y(1 : length(y) / 2 + 1); % フーリエ変換結果の後半分を破棄
f = linspace(0, 22050, length(Y)); % 周波数軸を作成
figure(2);
subplot(1,4,1);
plot(f, abs(Y)); % フーリエ変換結果を図示
xlim([0,2000]);
title('FOURIER');

% 短時間フーリエ変換
y1 = y(1 * fs + (0 : 4095)); % 開始から1 秒後の4096 点を取り出す
Y1 = fft(y1)/length(y1); % フーリエ変換
Y1 = Y1(1 : length(y1) / 2 + 1); % フーリエ変換結果の後半分を破棄
f = linspace(0, 22050, length(Y1)); % 周波数軸を作成
subplot(1,4,2);
plot(f, abs(Y1)); % フーリエ変換結果を図示
xlim([0,2000]);
title('NO-WINDOW');

% 短時間フーリエ変換

%(HANNING WINDOW)
y1 = y(1 * fs + (0 : 4095)); % 開始から1 秒後の4096 点を取り出す
w = hann(4096); % ハニング窓作成
y1 = w .* y1; % 作成したハニング窓をかける
Y1 = fft(y1)/length(y1); % フーリエ変換
Y1 = Y1(1 : length(y1) / 2 + 1); % フーリエ変換結果の後半分を破棄
f = linspace(0, fs/2, length(Y1)); % 周波数軸を作成
subplot(1,4,3);
plot(f, abs(Y1)); % フーリエ変換結果を図示
xlim([0,2000]);
title('HANNING');

%(HAMMING WINDOW)
y1 = y(1 * fs + (0 : 4095)); % 開始から1 秒後の4096 点を取り出す
w = hamming(4096); % ハニング窓作成
y1 = w .* y1; % 作成したハニング窓をかける
Y1 = fft(y1)/length(y1); % フーリエ変換
Y1 = Y1(1 : length(y1) / 2 + 1); % フーリエ変換結果の後半分を破棄
f = linspace(0, fs/2, length(Y1)); % 周波数軸を作成
subplot(1,4,4);
plot(f, abs(Y1)); % フーリエ変換結果を図示
xlim([0,2000]);
title('HAMMING');

警告: インデックスとして用いるコロン(:)を使用した演算では、整数を指定する必要があります

SPECTROGRAMのハニング窓の違いによる時間分解能と周波数分解能の変化の様子

figure(3);
spectrogram(y,hann(256),128,256,fs,'yaxis');
title('SPECTROGRAM with HANNING(256)');

figure(4);
spectrogram(y,hann(1024),512,1024,fs,'yaxis');
title('SPECTROGRAM with HANNING(1024)');

Published with MATLAB® 7.3