%% SNRの計算
clear all;

% パラメータ設定・取得 (org.wav:原曲)
N = 65536;      % サンプル数
ofst = 20000;   % 開始点からのオフセット
[dummy,FS,NBITS,OPTS] = wavread('org.wav',0);

% 音楽データ読み込み (ステレオ想定、org.wav:原曲、mp3.wav:MP3復号曲)
[oWAV,FS,NBITS,OPTS] = wavread('org.wav',[ofst+1,ofst+N]);
[cWAV,FS,NBITS,OPTS] = wavread('mp3.wav',[ofst+1,ofst+N]);

% FFT＆SNR計算(Leftのみ)
foWAV = abs(fft(oWAV(:,1)));
fcWAV = abs(fft(cWAV(:,1)));
fdif = foWAV-fcWAV;
fpow = mean( foWAV.^2 );
fsigma = mean( fdif.^2 );
fsnr = -10.0 * log10( fsigma / fpow )

% 波形、FFT結果表示
subplot(3,1,1);
plot(1:N, oWAV(:,1), 1:N, cWAV(:,1));
subplot(3,1,2);
plot(1:N, foWAV, 1:N, fcWAV);
%plot(1:N, log10(foWAV), 1:N, log10(fcWAV));
subplot(3,1,3);
plot(1:N, fdif, 'r');
%plot(1:N, log10(abs(fdif)), 'r');


