clear all;
[y, Fs, NBITS] = wavread('Test.wav');
y = filter([1 -0.97],1,y);                      % プリエンファシス
Fs = Fs * 1/2;                                  % サンプリングレートを1/2倍に（ピッチ下降）

frameSize = 0.025;								% フレーム長：0.025秒（25ms）
frameShift = 0.010;								% フレームシフト長：0.010秒（10ms）
frameSizeSample = fix( Fs * frameSize );		% フレーム長：サンプル換算
frameShiftSample = fix( Fs * frameShift );		% フレームシフト長：サンプル換算
maxFrame = fix((length(y)-(frameSizeSample-frameShiftSample))/frameShiftSample)-1;
startFrame = 1;									% フレームの開始サンプル番号
endFrame = startFrame + frameSizeSample - 1;	% フレームの終了サンプル番号

% 再生速度を何倍にするか（正の整数で指定して下さい）
speechRate = 2;

% 修正後の音声データの格納用
y2 = zeros(frameSizeSample, 1);
for countFrame = 1 : 1 : maxFrame
    thisData = y(startFrame : endFrame);
    window = hanning(frameSizeSample);  		% ハニング窓を生成
    thisData = thisData .* window;      		% 切り出したフレームにハニング窓をかける
    
    % countFrame が speechRate で割り切れる数の時だけ、切り出したフレームを【足し合わせる】
    if mod(countFrame, speechRate) == 0
        % フレームシフト部の足しあわせ
    	y2_shift = y2( length(y2)-frameShiftSample+1 : length(y2)) + thisData(1 : frameShiftSample);
        % フレームの追加
    	y2 = [y2(1:length(y2)-frameShiftSample) ; y2_shift ; thisData( frameShiftSample+1 : length(thisData) )];
    end

    startFrame = startFrame + frameShiftSample;
    endFrame = startFrame + frameSizeSample - 1;
end

sound(y2, Fs);                             			% 再生
wavwrite(y2, Fs, 'handel_pitchUp.wav');

%LSEE-MSTFTMアルゴリズム
%【各フレームの振幅スペクトル】の計算
startFrame = 1; endFrame = startFrame + frameSizeSample - 1;
% 音声の長さが変わったので最大フレーム数も計算し直す
maxFrame = fix((length(y2)-(frameSizeSample-frameShiftSample))/frameShiftSample)-1;
ModifiedSTFTM = zeros(frameSizeSample*2, maxFrame); % 各フレームの振幅スペクトルの格納用
for countFrame = 1 : 1 : maxFrame
    thisData = y2(startFrame : endFrame);
    window = hanning(frameSizeSample); thisData = thisData .* window;
    fftsize = frameSizeSample * 2;                  % FFT次数
    dft = fft(thisData, fftsize);                   % フーリエ変換
    Adft = abs(dft);                                % 振幅スペクトル
    ModifiedSTFTM(:, countFrame) = Adft;            % 各フレームの振幅スペクトルを記録
    startFrame = startFrame + frameShiftSample;
    endFrame = startFrame + frameSizeSample - 1;
end
%【各フレームの位相スペクトル】の計算
wgn_sf = Fs; wgn_n = length(y2); rand('state', sum(100 * clock));
y3 = randn(1, wgn_n)';                              % 初期信号（ホワイトノイズ）の生成
% 50回のループ
for count_i = 1 : 50
    y_new = zeros(length(y3), 1);
    % 各フレームごとに解析
    startFrame = 1; endFrame = startFrame + frameSizeSample - 1;
    for count_frame = 1 : maxFrame
        thisData = y3(startFrame : endFrame);
        window = hanning(frameSizeSample);
        thisData = thisData .* window;
        dft_LSEEMSTFTM = fft(thisData, fftsize);
        %【data_LSEEMSTFTM の位相スペクトル】
        Angdft_LSEEMSTFTM = atan2(imag(dft_LSEEMSTFTM), real(dft_LSEEMSTFTM));
        %【振幅スペクトルとdata_LSEEMSTFTM の位相スペクトルから音声波形を作る】
        sound_fft_LSEEMSTFTM = ModifiedSTFTM(:,count_frame) .* exp(Angdft_LSEEMSTFTM*j);
        % 逆フーリエ変換
        sound_part_LSEEMSTFTM = real(ifft(sound_fft_LSEEMSTFTM, fftsize));
        sound_part_LSEEMSTFTM = sound_part_LSEEMSTFTM(1:fftsize/2);
        %音声を足しあわせ
        y_new(startFrame:endFrame) = y_new(startFrame:endFrame) + sound_part_LSEEMSTFTM;
        startFrame = startFrame + frameShiftSample;
        endFrame = startFrame + frameSizeSample - 1;
    end
    y3 = y_new;
end

sound(y3, Fs);                                      % 再生（LSEE_MSTFTMアルゴリズム適用後）
wavwrite(y3, Fs, 'Test_pitchDOWN_sampleratedown.wav');