addpath(genpath(\'A:\SpeechIntelligibility\matlab\matlab_system\'));
clear
close all
fs = 16000;
N = 128;
N_frame = N;
f0 = fs/N * (N/2)
f1 = fs/N * (N/4-2);
f2 = fs/N * (N/4-3);
l = ceil(2*fs/N);
% l = 1;
t = 0:1/fs:(l*N-1)/fs;
% data = cos(2*pi*f0*t)+cos(2*pi*f1*t)+cos(2*pi*f2*t);
data = cos(2*pi*f0*t);
% data = randn(length(t),1);
data = data(:);
% w1 = fft(data);
% figure(\'position\',[0 0 400 275]); plot(abs(w1),\'-o\')
N_fft = N_frame;
N_shift = N_frame/4;
win1 = hann(N_frame,\'periodic\');
% win1 = rectwin(N_frame);
win = WinSqrt(win1,N_shift);
[a1, X] = stft(data, N_frame, N_fft, N_shift, win);
w1 = inverse_stft(X, N_frame, N_fft, N_shift, win);
% spclab(fs, a1, w1);
% figure(\'position\',[0 0 400 275]); plot(mean(abs(X(1:end/2+1,:)),2),\'-o\')
%% Anti-aliasing
Y = X;
Y_half = Y(1:end/2+1,:);
p = 2;
y1 = ifft(X);
y1 = lp_time(y1, (N_frame/2)/p);
%% Antialiasing realized by sinc function.
% w = (0:N/2)*0.5;
% sinc1 = N/2*sinc(w);
% sinc1 = sinc1(:);
% % win2 = hann(N+1,\'symmetric\');
% win2 = gausswin(N+1,3.5);
% win2_half = win2(N/2+1:end);
% H1_half = sinc1 .* win2_half;
% H1 = spec_half2full(H1_half);
% h1 = ifft(H1);
% % figure; plot(h1)
% h1_array = repmat(h1, 1, size(y1,2));
% y1 = y1 .* h1_array;
X = fft(y1);
%%
X_half = X(1:N_fft/2+1,:);
% X1_half = zeros(size(X_half));
% % X1_half(end-5-10:end-10,:) = X_half(end-5:end,:);
% shift = 30;
% X1_half(1+shift:20+shift,:) = X_half(1:20,:);
% X1_half = zeros(size(X_half));
% % tmp1 = X_half(2,1);
% % X1_half(50,:) = ones(1,size(X1_half,2))*tmp1;
% X1_half(2,:) = X_half(64,:);
X1_half = X_half;
% X1_half(2,:) = 0;
% X1_half = X_half(1:p:end,:);
% X1_half(end+1:N_fft/2+1,:) = 0;
%% DEBUG
figure(\'position\',[0 0 400 275]); plot(mean(abs(X_half),2),\'-o\')
figure(\'position\',[0 0 400 275]); plot(mean(abs(X1_half),2),\'-*\')
%% Use original phase
% X1_half_abs = abs(X1_half);
% X_half_phase = angle(Y_half);
% X1_half = X1_half_abs .* exp( 1i*X_half_phase );
%%
X1 = spec_half2full(X1_half);
x1 = inverse_stft(X1, N_frame, N_fft, N_shift, win);
%% TEST
x2 = ifft(X1);
%%
% sound(x1,fs);
% Method 1
% p_a1 = pwelch(w1, rectwin(N));
% p_x1 = pwelch(x1, rectwin(N));
% Method 2
% p_a1 = pwelch(w1);
% p_x1 = pwelch(x1);
% figure; hold on; plot(10*log10(p_a1)); plot(10*log10(p_x1)); hold off
% Method 3
win2 = hann(N_frame,\'periodic\');
[p_a1,faxis1] = pwelch(w1,win,N_frame/8,N_frame,fs);
[p_x1, faxis2] = pwelch(x1,win,N_frame/8,N_frame,fs);
figure; hold on; plot(faxis1, 10*log10(p_a1+eps)); plot(faxis2,10*log10(p_x1+eps)); hold off
figure; plot(x1)
spclab(fs, a1, x1);