% 参考代码： Matlab 信号处理的基本操作，生成信号，从时域和频域观察信号，对信号滤波
% 每个功能演示代码块都使用了 if( ) ... end 来封装。通过使能不同位置的 if()判断变量，启用不同的代码块。
% 对于不熟悉的Matlab函数，在Matlab命令行输入 doc xxx 即可阅读帮助文档，例如 doc fir2  
% 作者：中国传媒大学 杜伟韬 wesignal@qq.com
% 本代码在Matlab2020b下测试通过
% 欢迎中传在校生同学报名参加当年的校内电子设计竞赛进而参加全国大赛
% https://ecdav.cuc.edu.cn/37/list.htm


clear;clc;close all;
% 生成 正弦波信号并且观察波形和FFT频谱 时频光谱图
if(1)
    fs = 8E3 ;         % samplerate in Hz
    f0 = 1500;         % signal frequency in Hz
    N =  200;         % Num of Points
    idx = [0:N-1];   % index vector
    t   =  idx/fs; % time vector
    x = sin(2*pi*f0*t);
    figure;plot(t,x);grid on;title(['sine signal, f0:' num2str(f0) ' Hz,  fs:', num2str(fs), ' Hz']);
    f = idx / N * fs;
    y = abs(fft(x));
    y_dB = 20*log10(y);
    figure;plot(f, y_dB);grid on; title(['Sine signal ' num2str(N) ' points fft, in dB scale, fs:' num2str(fs), ' Hz'] );
    spg_win_len = 128; spg_overlap_len = 100; spg_fft_len = 128;
    figure;spectrogram(x, spg_win_len,spg_overlap_len,spg_fft_len,fs); title('Sine Signal Spectrogram');
end
% 生成 多音（Multi-Tone）正弦波信号，观察波形和FFT频谱 时频光谱图
if(0)
    fs = 8E3 ;         % samplerate in Hz
    f0 =  300;         % signal frequency in Hz
    f1 =  500;         % signal frequency in Hz
    f2 =  800;         % signal frequency in Hz
    f3 = 1500;         % signal frequency in Hz
    N =  800;          % Num of Points
    f_vec = [f0 f1 f2 f3 ];
    idx = [0:N-1];     % index vector
    t   =  idx/fs;     % time vector
    x0 = sin(2*pi*f0*t);
    x1 = sin(2*pi*f1*t);
    x2 = sin(2*pi*f2*t);
    x3 = sin(2*pi*f3*t);
    x = (x0+x1+x2+x3)/4;
    figure;plot(t,x);grid on;title([num2str(f_vec) 'Hz,  multi-tone sine wave, fs:', num2str(fs), ' Hz']);
    f = idx / N * fs;
    y = abs(fft(x));
    y_dB = 20*log10(y);
    figure;plot(f, y_dB);grid on; title(['Multi-tone signal' num2str(N) ' points fft, in dB scale, fs:' num2str(fs), ' Hz'] );
    spg_win_len = 128; spg_overlap_len = 100; spg_fft_len = 128;
    figure;spectrogram(x, spg_win_len,spg_overlap_len,spg_fft_len,fs); title('Multi-tone Signal Spectrogram');
end


% 生成一个频率阶跃的正弦波信号并观察波形、FFT频谱和时频光谱图

if(0)
    fs = 8E3 ;         % samplerate in Hz
    f0 =  150;         % signal frequency in Hz
    f1 = 1500;         % signal frequency in Hz
    N0 =  400;         % Num of Points
    N1 =  400;         % Num of Points
    N = N0+N1;
    f_vec = [f0 f1];
    idx0= [0:N0-1];   % index vector
    idx1= [0:N1-1];   % index vector
    idx = [0: N-1];   % index vector
    t0  =  idx0/fs;   % time vector
    t1  =  idx1/fs;   % time vector
    t   =  idx/fs ;   % time vector
    x0= sin(2*pi*f0*t0);
    x1= sin(2*pi*f1*t1);
    x = [x0 x1];
    figure;plot(t,x);grid on;title([num2str(f_vec) 'Hz,  freq step sine wave, fs:', num2str(fs), ' Hz']);
    f = idx / N * fs;
    y = abs(fft(x));
    y_dB = 20*log10(y);
    figure;plot(f, y_dB);grid on; title(['Freq step Sine signal ' num2str(N) ' points fft, in dB scale, fs:' num2str(fs), ' Hz'] );
    spg_win_len = 128; spg_overlap_len = 100; spg_fft_len = 128;
    figure;spectrogram(x, spg_win_len,spg_overlap_len,spg_fft_len,fs); title('Freq step sine Signal Spectrogram');
end


% 设计FIR滤波器观察频响
if (0)
    fs    = 8E3 ;   % samplerate in Hz
    f_cut = 1E3 ;   % filter cut-off freq in Hz 
    n_odr = 64  ;   % fir order   
    f_set = [0 f_cut f_cut  fs/2]/(fs/2); % freq vec for fir
    m_set = [1 1 0 0];
    b_fir = fir2(n_odr,f_set,m_set);

    freqz(b_fir,1,[],fs) ; title('FIR Freq response');
end

% 生成多音正弦，用FIR滤波器对多音正弦波信号滤波，观察时域波形和FFT频谱
if(0)
    fs = 8E3 ;         % samplerate in Hz
    f0 =  300;         % signal frequency in Hz
    f1 =  500;         % signal frequency in Hz
    f2 =  800;         % signal frequency in Hz
    f3 = 1500;         % signal frequency in Hz
    N =  800;          % Num of Points
    f_vec = [f0 f1 f2 f3 ];
    idx = [0:N-1];     % index vector
    t   =  idx/fs;     % time vector
    x0 = sin(2*pi*f0*t);
    x1 = sin(2*pi*f1*t);
    x2 = sin(2*pi*f2*t);
    x3 = sin(2*pi*f3*t);
    x = (x0+x1+x2+x3)/4;
    figure;plot(t,x);grid on;title([num2str(f_vec) 'Hz,  multi-tone sine wave, fs:', num2str(fs), ' Hz']);
    f = idx / N * fs;
    y = abs(fft(x));
    y_dB = 20*log10(y);
    figure;plot(f, y_dB);grid on; title(['Multi-tone signal' num2str(N) ' points fft, in dB scale, fs:' num2str(fs), ' Hz'] );
    spg_win_len = 128; spg_overlap_len = 100; spg_fft_len = 128;
    figure;spectrogram(x, spg_win_len,spg_overlap_len,spg_fft_len,fs); title('Multi-tone Signal Spectrogram');
    
    f_cut = 600 ;   % filter cut-off freq in Hz 
    n_odr = 64  ;   % fir order   
    f_set = [0 f_cut f_cut  fs/2]/(fs/2); % freq vec for fir
    m_set = [1 1 0 0];
    b_fir = fir2(n_odr,f_set,m_set);
    figure; freqz(b_fir,1,[],fs);title('FIR Freq response');
    
    out = filter(b_fir,1,x); 
    figure;plot(t,out);grid on;title([' FILTERED, ' num2str(f_vec) 'Hz, multi-tone sine wave, fs:', num2str(fs), ' Hz']);
    yout = abs(fft(out));
    yout_dB = 20*log10(yout);
    figure;plot(f, yout_dB);grid on; title([' FILTERED, Multi-tone signal' num2str(N) ' points fft, in dB scale, fs:' num2str(fs), ' Hz'] );
    figure;spectrogram(out, spg_win_len,spg_overlap_len,spg_fft_len,fs); title('FILTERED, Multi-tone Signal Spectrogram');

end

% 用FIR滤波器对频率阶跃的正弦波信号滤波，观察时域波形和FFT频谱
if(0)
    fs = 8E3 ;         % samplerate in Hz
    f0 =  150;         % signal frequency in Hz
    f1 = 1500;         % signal frequency in Hz
    N0 =  400;         % Num of Points
    N1 =  400;         % Num of Points
    N = N0+N1;
    f_vec = [f0 f1];
    idx0= [0:N0-1];   % index vector
    idx1= [0:N1-1];   % index vector
    idx = [0: N-1];   % index vector
    t0  =  idx0/fs;   % time vector
    t1  =  idx1/fs;   % time vector
    t   =  idx/fs ;   % time vector
    x0= sin(2*pi*f0*t0);
    x1= sin(2*pi*f1*t1);
    x = [x0 x1];
    figure;plot(t,x);grid on;title([num2str(f_vec) 'Hz,  freq step sine wave, fs:', num2str(fs), ' Hz']);
    f = idx / N * fs;
    y = abs(fft(x));
    y_dB = 20*log10(y);
    figure;plot(f, y_dB);grid on; title(['Freq step Sine signal ' num2str(N) ' points fft, in dB scale, fs:' num2str(fs), ' Hz'] );
    spg_win_len = 128; spg_overlap_len = 100; spg_fft_len = 128;
    figure;spectrogram(x, spg_win_len,spg_overlap_len,spg_fft_len,fs); title('Freq step sine Signal Spectrogram');
    
    f_cut = 600 ;   % filter cut-off freq in Hz 
    n_odr = 64  ;   % fir order   
    f_set = [0 f_cut f_cut  fs/2]/(fs/2); % freq vec for fir
    m_set = [1 1 0 0];
    b_fir = fir2(n_odr,f_set,m_set);
    figure; freqz(b_fir,1,[],fs);title('FIR Freq response');
    
    out = filter(b_fir,1,x); 
    figure;plot(t,out);grid on;title([' FILTERED, ' num2str(f_vec) 'Hz, Freq step sine wave, fs:', num2str(fs), ' Hz']);
    yout = abs(fft(out));
    yout_dB = 20*log10(yout);
    figure;plot(f, yout_dB);grid on; title([' FILTERED, Freq step sine signal' num2str(N) ' points fft, in dB scale, fs:' num2str(fs), ' Hz'] );
    figure;spectrogram(out, spg_win_len,spg_overlap_len,spg_fft_len,fs); title('FILTERED, Freq step sine Signal Spectrogram');

end