clear all

% channel 1: target point
% channel 2: left ear
% channel 3: right ear

fn_data = 'LeftHandel_Right40kHz_Record.wav';
[sig, fs] = audioread(fn_data);

% ch = 1;
% sig.y = sig.y(:, ch);
for ch = 1:3
    y_now = sig(:, ch);
    y_now = highpass(lowpass(y_now, 0.01), 50, fs);
    % p = audioplayer(y_now*1e3, fs); play(p);

    %% 
    figure;
    pwelch(y_now, hanning(fs), floor(fs/2), fs, fs);
    xlim([0,4])
    
    %%
    win_len = 4096; % 帧长
    win = hamming(win_len,'periodic');% 窗口函数（汉明窗）
    overlap = win_len/2; %帧移（一般为帧长的一半）
    fft_len = win_len; %做DFT的点数，一般和帧长一样
    t = (1/fs) * (0:numel(y_now)-1);%波形图的时间（单位为秒）
    
    sig_ampl = 2e3;
    t_shift = 0.4;
    
    %% 
    y = y_now * sig_ampl;
    figure;
    plot(t, y); 
    title('Time signal');
    xlim([0,8.9])
    ylim([-1,1])
    
    %%
    figure;
    [~, f, t, pd] = spectrogram(y, win, overlap, fft_len, fs,'yaxis');
    pcolor(t-t_shift, f, 10*log10(abs(pd)))
    clb = colorbar;
    colormap(JiaxinColor('plasma'))
    shading interp
    set(gca, 'fontsize', 20)
    xlabel('Time (s)');
    ylabel('Frequency (Hz)')
    pbaspect([4,1,1])
    ylim([0,4e3])
    xlim([0,8.9])
    set(gca, 'xtick', [0:1:9])
    set(gca, 'ytick', [0:2e3:4e3]);
    set(clb, 'ytick', [-100, -75, -50]);
    clim([0, 50] - 100);
    set(gcf, 'position', [248.2	267	926.8	420])
end