%% 產生FM信號並繪製時域波形
fm = 500; % 調製信號頻率
Am = 0.5; % 調製信號幅度
fc = 5e3; % 載波頻率
Ac = 1; % 載波幅度
kf = 10000*pi; % 調頻靈敏度
fs = 75e3; % 採樣率
N = 3000; % 樣點總數
t = (0:N-1)'/fs; % 時間t
m_t = Am*sin(2*pi*fm*t); % 調製信號
phi_t = kf*cumsum(m_t)/fs; % 相位積分
s_t = cos(2*pi*fc*t + phi_t); % 已調信號
plot(t, s_t , 'b'); % 繪波形
xlabel('time');
xlim('auto');
ylabel('amplitude');
ylim('auto');

%% 繪製頻譜圖
L = length(s_t); % 取得序列長度 
u = fftshift(fft(s_t )); % 離散傅里葉變換，求頻譜
w = (0:L-1)'*fs/L - 1/2*fs; % 橫座標-頻率
figure(2)
plot(w, abs(u)); % 繪製頻譜圖 
grid on; 
xlabel('frequency(Hz)'); 
ylabel('magnitude(dB)');

%% 貝塞爾函數值
b=500*N/fs;%計算頻域上間隔500Hz對應的點數
c=1500+10*b+1;%計算fc在頻域上對應的位置
for i=0:1:10
   a(i+1)=abs(u(c+b*i)/N*2);%計算fc+i*fm對應的頻譜值
   bessel(i+1)=besselj(i,5);%計算對應的貝塞爾函數值，mf=5
end
a
bessel

%FM信號非相干解調
fc = 5e3; % 載波頻率
fs = 100e3; % 採樣率
N = length(s_t); % 樣點數
[lpf_b,lpf_a] = butter(10, (fc/2)/(fs/2)); % 設計低通濾波器
t = (0:N-1)'/fs; % 時間t
r_t = s_t;
r_d_t = [0;diff(r_t)]; % 求微分
r_e_t = abs(hilbert(r_d_t)); % 包絡檢波
demod_t = filter(lpf_b, lpf_a, r_e_t); % 濾波
plot(t, demod_t , 'b'); % 繪圖

%% 調頻立體聲廣播解調 
clear all;clc;close all
load 'fm_cap.mat';  %載入FM信號
fc = 200e3; % 載波頻率
fs = 2e6; % 採樣率
r_t = fm_cap;
N = length(r_t); % 樣點數
[lpf_b,lpf_a] = butter(10, (fc/2)/(fs/2)); % 設計低通濾波器
t = (0:N-1)'/fs; % 時間t
%% 解調
r_d_t = [0;diff(r_t)]; % 求微分
r_e_t = abs(hilbert(r_d_t)); % 包絡檢波
demod_t = filter(lpf_b, lpf_a, r_e_t); % 濾波
plot(t, demod_t , 'b'); % 繪圖
%% 信號抽取，使聲卡能夠播放音頻
demod_t2 = resample(demod_t,48e3,2e6);%按比例降低採樣率
demod_t2=demod_t2/max(demod_t2);%歸一化，限制音頻幅度範圍，否則無法播放
fs2=48e3;%新採樣率
N2 = length(demod_t2); % 新樣點數
t2 = (0:N2-1)'/fs2; % 新時間t
plot(t2, demod_t2 , 'b'); % 繪圖
sound(demod_t2,fs2);%播放音頻