Matlabでオーディオファイルを扱う

パスを指定

Macの場合

cd('/Users/＊＊＊/＊＊＊/＊＊＊')

プログラムの最後に、Mファイルを保存しているディレクオリに戻る記述をしておくと便利。

オーディオの読み込み

s = audioread('sample.wav');
[s,fs] = audioread('sample.wav');

オーディオの書き出し

名前を指定してオーディオをwav形式で書き出します

filename = 'name.wav';
audiowrite(filename,y,Fs);
clear y Fs

スペクトログラムを定義する関数

function [S, time, frequency] = Spectrogram_(s, fs, window_size, shift_size)
dft_size = window_size;
length_of_s = length(s);
number_of_frame = floor((length_of_s -(window_size - shift_size))/shift_size);
x = zeros(1, window_size);
w = HanningWindow_(window_size);
S = zeros(dft_size/2+1, number_of_frame);
for frame = 1:number_of_frame
    offset = shift_size*(frame - 1);
    for n = 1:window_size
        x(n) = s(offset+n)*w(n);
    end
    X = fft(x, dft_size);
    for k = 1:dft_size/2+1
        S(k, frame) = 20*log10(abs(x(k)));
    end
end
time = zeros(1, number_of_frame);
for frame = 1:number_of_frame
    time(frame) = (frame - 1)*shift_size/fs;
end
frequency = zeros(1, dft_size/2+1);
for k = 1:dft_size/2+1
    frequency(k) = (k-1)*fs/dft_size;
end

DFTを定義する関数

離散フーリエ変換の定義式

function [X_real,X_imag]=DFT_(x_real,x_imag,dft_size)
X_real=zeros(1,dft_size);
X_imag=zeros(1,dft_size);

for k=1:dft_size
    for n=1:dft_size
    w_real=cos(2*pi*(k-1)*(n-1)/dft_size);
    w_imag=-sin(2*pi*(k-1)*(n-1)/dft_size);
    X_real(k)=X_real(k)+w_real*x_real(n)-w_imag*x_imag(n);
    X_imag(k)=X_imag(k)+w_real*x_imag(n)-w_imag*x_real(n);
    end
end

IDFTを定義する関数

逆離散フーリエ変換の定義式

function [x_real,x_imag]=IDFT_(X_real,X_imag,dft_size)
x_real=zeros(1,dft_size);
x_imag=zeros(1,dft_size);

for n=1:dft_size,
    for k=a:dft_size,
    w_real=cos(2*pi*(k-1)*(n-1)/dft_size);
    w_imag=sin(2*pi*(k-1)*(n-1)/dft_size);
    x_real(n)=x_real(n)+w_real*X_real(k)-w_imag*X_imag(k);
    x_imag(n)=x_imag(n)+w_real*X_imag(k)-w_imag*X_real(k);
    end
    x_real(n)=x_real(n)/dft_size;
    x_imag(n)=x_imag(n)/dft_size;
end

FFTを定義する関数

function [X_real,X_imag]=FFT_(x_real,x_imag,dft_size)
X_real=zeros(1,dft_size);
X_imag=zeros(1,dft_size);
number_of_stage=log2(dft_size);
for stage=1:number_of_stage
    for i=1:power(2,(stage-1))
        for j=1:power(2,(number_of_stage-stage))
            n=power(2,(number_of_stage-stage+1))*(i-1)+j;
            m=power(2,(number_of_stage-stage))+n;
            r=power(2,(stage-1))*(j-1);
            a_real=x_real(n);
            a_imag=x_imag(n);
            b_real=x_real(m);
            b_imag=x_imag(m);
            c_real=cos((2*pi*r)/dft_size);
            c_imag=-sin((2*pi*r)/dft_size);
            if stage<number_of_stage
                x_real(n)=a_real+b_real;
                x_imag(n)=a_imag+b_imag;
                x_real(m)=(a_real-b_real)*c_real-(a_imag-b_imag)*c_imag;
                x_imag(m)=(a_imag-b_imag)*c_real+(a_real-b_real)*c_imag;
            else
                x_real(n)=a_real+b_real;
                x_imag(n)=a_imag+b_imag;
                x_real(m)=a_real-b_real;
                x_imag(m)=a_imag-b_imag;

            end
        end
    end
end
number_of_bit=log2(dft_size);
for k=1:dft_size
    index=dec2bin((k-1),number_of_bit);
    for position_of_bit=1:floor(number_of_bit/2)
        bit=index(position_of_bit);
        index(position_of_bit)=index(number_of_bit-position_of_bit+1);
        index(number_of_bit-position_of_bit+1)=bit;
    end
    X_real(bin2dec(index)+1)=x_real(k);
    X_imag(bin2dec(index)+1)=x_imag(k);
end

IFFTを定義する関数