% 参考代码：FM波段信号 滤波、解调算法
% 作者：中国传媒大学 杜伟韬 wesignal@qq.com
% 本代码在Matlab2020b下测试通过
% 数据文件来自 http://www.dudulab.net/radio_FM/content.htm
% 欢迎中传在校生同学报名参加当年的校内电子设计竞赛进而参加全国大赛
% https://ecdav.cuc.edu.cn/37/list.htm


function main_fm_dec()
% CF 和 fs 参数，请设定为 网页频谱图中给出的参数
  CF = 103.9E6; % 中心频率 Hz
  fs = 6.144E6;  % Hz
  nbits = 16;   % 比特
  fname_rd_iq = 'fm1.dat';
  max_t_len = 5;    % 秒
  iq_rd = open_rd_iq_file( fs, nbits, max_t_len, fname_rd_iq);

  in = iq_rd;
  % design halfband filter coeff
  b_hbf = firhalfband('minorder',.45,0.0001);
  % generate hbf object
  firdecim = dsp.FIRDecimator(2,b_hbf);
  %  filter the signal
  disp('DECIM 2X');
  in_dcm2 = firdecim(in);
  disp('DECIM 4X');
  in_dcm4 = firdecim(in_dcm2);
  disp('DECIM 8X');
  in_dcm8 = firdecim(in_dcm4);
  disp('DECIM 16X');
  in_dcm16 = firdecim(in_dcm8);
  disp('DECIM 32X');
  in_dcm32 = firdecim(in_dcm16);
 
  in_base = in_dcm32; 
  fs_b    = fs/32;
  
  disp('GET POWER SPECTRUM');
  % use: [pxx,f] = pwelch(x,window,noverlap,f,fs)
  figure;pwelch(in,8192,[],[],fs,'centered','power');title(['IF Sampled FM signal, sample rate :', num2str(fs/1E3),' KHz, Center Freq : ', num2str(CF/1E6), ' MHz']);grid on;
  figure;pwelch(in_base,8192,[],[],fs_b,'centered','power'); title(['Base Band FM signal, sample rate :', num2str(fs_b/1E3), ' KHz, Center Freq : ', num2str(CF/1E6), ' MHz']); grid on;

  disp('MONO FM DEC ');
  vec_ang = angle(in_base);
  vec_a_df = diff(unwrap(vec_ang))/pi; % must use unwrap function, in order to make the phase difference smaller than pi
  vec_mpx = [0;vec_a_df];
  
  figure;pwelch(vec_mpx,8192,[],[],fs_b,'onesided' ,'power');title(['FM Demod signal, sample rate :', num2str(fs_b/1E3), ' KHz']);

  
  b_au_hbf = firhalfband('minorder',.4,0.0001);
  fir_au_decim = dsp.FIRDecimator(2,b_au_hbf);
  
  disp('BASE BAND DECIM '); 
  vec_mpx_d2 = fir_au_decim(vec_mpx);
  vec_mpx_d4 = fir_au_decim(vec_mpx_d2);
  au_out = vec_mpx_d4;
  fs_au = fs_b/4;
  disp('AUDIO WRITE OUT'); 
  audiowrite('au_out.wav',au_out,fs_au,'BitsPerSample',16);


  disp('main, DONE');

end

function iq_out = open_rd_iq_file( fs, nbits, max_t_len, fname);

    max_rd_sample =  fs*max_t_len*2;
    max_rd_byte = fix(max_rd_sample*(nbits/8));    
    % 测试文件大小
    fid = fopen(fname);
    if(fid == -1)
        disp('OPEN FILE failed');        
    end
    fseek(fid,0,'eof'); % 到文件结尾
    fsize_byte = ftell(fid); %fsize是文件大小
    rd_byte = max_rd_byte; 
    if(fsize_byte < max_rd_byte)
        rd_byte = fsize_byte;
    end
    rd_sample = fix(rd_byte*8/nbits);

    fseek(fid,0,'bof'); % 到文件开头
    v_rd = fread(fid,rd_sample,'int16');
    v_iq = v_rd(1:2:rd_sample-1) + j*v_rd(2:2:rd_sample);

    fclose(fid);
    iq_out = v_iq;
    if(0)
      disp('INTO DEBUG LOOP open_rd_iq_file()');
      while(1)
      end
    end


end