前前后后花了三天时间
其实第一天就把fmmod 和 fmdemod做出来了
只是感觉很多东西 跟我理解的匹配不上
一个是带宽 一个是最大频偏
最后实现了fmmod 和pmmod实现FM调制 再使用fmdemod解调
关于带宽 现在还不是很理解 因为画出来的频谱非常非常集中 范围并不大 为什么带宽需要这么大（做了一个更小带宽（能量占比也达到了0.99）的带通 处理调制信号 再解调 效果是差一些 ）
关于瞬时频率 instfreq不能用到这个地方 因为不是单分量 至于为什么 还是不知道
做出来fmmod 和pmmod的信噪比 有一些差距 不知道是不是正常的
%Principle of Communications Topic 1 FM %Contributed by LiuYH clc;clear;close all; %读取音频文件 [original_signal, fs] = audioread('testaudio.wav'); %幅度归一化处理 max_amplitude = max(abs(original_signal)); original_signal = original_signal / max_amplitude; %设置FM调制参数 B = 8000; %基带信号带宽 计算见BW.m Kfm = 2000; %频偏常数 fc = 11000; %载波频率 beta = Kfm/B; %调频指数/频率调制指数 %时间向量 绘图用 t = linspace(0, length(original_signal)/fs, length(original_signal)); len = round(length(original_signal)/300); %只画一部分长度 % 画出original_signal的时域图 figure; plot(t, original_signal); xlabel('Time (s)'); ylabel('Amplitude'); title('original signal 时域图'); %画出original_signal的频谱图 N=length(original_signal); figure; f = (0:N-1)*fs/N - fs/2; plot(f,20*log10(abs(fftshift(fft(original_signal))))); xlabel('Frequency (Hz)'); ylabel('Amplitude (dB)'); title('original signal 频域图'); %验证基带信号带宽B passband = [0 7880]; % 通带频率范围 % 使用傅里叶变换将modulated_signal转换为频域表示 与前面画频谱的时候是一样的操作 但是不在乎复杂度 hhh original_signal_spectrum = fftshift(fft(original_signal)); % 确定要保留的频率范围，将其他频率分量置零 f1 = (-fs/2 : fs/length(original_signal) : fs/2 - fs/length(original_signal)); original_signal_spectrum(abs(f1) < passband(1) | abs(f1) > passband(2)) = 0; filtered_signal_spectrum = original_signal_spectrum; % 使用逆傅里叶变换将修改后的频域信号转换回时域表示 filtered_original_signal = ifft(ifftshift(filtered_signal_spectrum)); % 计算滤波前和滤波后信号的功率值 original_signal_energy = sum(abs(original_signal).^2); filtered_original_signal_energy = sum(abs(filtered_original_signal).^2); %将filtered_original_signal生成音频 % audiowrite('filtered_original_signal.wav',filt_signal,fs); %计算滤波前后的能量 fprintf('滤波前的能量：%f\n', original_signal_energy); fprintf('滤波后的能量：%f\n', filtered_original_signal_energy); fprintf('比值：%f\n', filtered_original_signal_energy/original_signal_energy); fprintf('\n'); % 绘制滤波后的信号频谱图 figure; plot(f,20*log10(abs(fftshift(fft(filtered_original_signal))))); title('filtered original signal 频域图'); xlabel('Frequency (Hz)'); ylabel('Amplitude (dB)'); %FM调制 fmmodulated_signal = fmmod(original_signal, fc, fs, Kfm); %使用PM调制实现FM调制 int_x = cumtrapz(t,original_signal); pmmodulated_signal = pmmod(2*pi*Kfm*int_x,fc,fs,1); %比较两种方法得到的调制信号 以FM调制的作为基准 signal_power_mod = mean(abs(fmmodulated_signal).^2); noise_power_mod = mean(abs(pmmodulated_signal - fmmodulated_signal).^2); snr_db_pmfm = 10 * log10(signal_power_mod / noise_power_mod); figure; plot(t(1:len),[pmmodulated_signal(1:len) fmmodulated_signal(1:len)]); legend('pmmod signal','fmmod signal'); %验证参数之最大频偏(分析瞬时频率) %不可行的方法 对单分量信号才可以使用instfreq %figure; %instfreq(modulated_signal,fs); %可行的方法 z = hilbert(fmmodulated_signal); instfrq = fs/(2*pi)*diff(unwrap(angle(z))); % audiowrite('recovered_signal_instfreq.wav',instfrq/10000-1.1,fs); figure; plot(t(2:end),instfrq); %验证参数之时域 %画出fmmodulated signal时域图 figure; plot(t, fmmodulated_signal); xlabel('Time (s)'); ylabel('Amplitude'); title('FMmodulated signal 时域图'); %画出pmmodulated signal时域图 figure; plot(t, pmmodulated_signal); xlabel('Time (s)'); ylabel('Amplitude'); title('PMmodulated signal 时域图'); %验证参数之频域 %画出fmmodulated signal频域图 figure; plot(f,20*log10(abs(fftshift(fft(fmmodulated_signal))))); title('FMmodulated signal 频域图'); xlabel('Frequency (Hz)'); ylabel('Amplitude (dB)'); %画出pmmodulated signal频域图 figure; plot(f,20*log10(abs(fftshift(fft(pmmodulated_signal))))); title('PMmodulated signal 频域图'); xlabel('Frequency (Hz)'); ylabel('Amplitude (dB)'); %验证参数之带宽 卡森公式 P_low = fc-B-B*beta; P_high = fc+B+B*beta; passband = [P_low P_high]; % 通带频率范围 % 使用傅里叶变换将fmmodulated_signal转换为频域表示 fmmodulated_spectrum = fftshift(fft(fmmodulated_signal)); pmmodulated_spectrum = fftshift(fft(pmmodulated_signal)); % 确定要保留的频率范围，将其他频率分量置零 f1 = (-fs/2 : fs/length(fmmodulated_signal) : fs/2 - fs/length(fmmodulated_signal)); fmmodulated_spectrum(abs(f1) < passband(1) | abs(f1) > passband(2)) = 0; filtered_fmmodulated_spectrum = fmmodulated_spectrum; pmmodulated_spectrum(abs(f1) < passband(1) | abs(f1) > passband(2)) = 0; filtered_pmmodulated_spectrum = pmmodulated_spectrum; % 使用逆傅里叶变换将修改后的频域信号转换回时域表示 filtered_fmmodulated_signal = ifft(filtered_fmmodulated_spectrum); filtered_pmmodulated_signal = ifft(ifftshift(filtered_pmmodulated_spectrum)); % 计算滤波前和滤波后信号的功率值 fmmodulated_energy = sum(abs(fmmodulated_signal).^2); pmmodulated_energy = sum(abs(pmmodulated_signal).^2); filtered_fm_energy = sum(abs(filtered_fmmodulated_signal).^2); filtered_pm_energy = sum(abs(filtered_pmmodulated_signal).^2); % 比较滤波前后的功能量 fprintf('滤波前的能量：%f\n', fmmodulated_energy); fprintf('滤波后的fmmod能量：%f\n', filtered_fm_energy); fprintf('比值为：%f\n',filtered_fm_energy/fmmodulated_energy); fprintf('滤波前的能量：%f\n', pmmodulated_energy); fprintf('滤波后的pmmod能量：%f\n', filtered_pm_energy); fprintf('比值为：%f\n',filtered_pm_energy/fmmodulated_energy); % 绘制滤波后的信号频谱图 figure; plot(f,20*log10(abs(fftshift(fft(filtered_fmmodulated_signal))))); xlabel('Frequency (Hz)'); ylabel('Amplitude (dB)'); title('filtered fmmod signal 频域图'); figure; plot(f,20*log10(abs(fftshift(fft(filtered_pmmodulated_signal))))); xlabel('Frequency (Hz)'); ylabel('Amplitude (dB)'); title('filtered pmmod signal 频域图'); %解调 recoverfm = fmdemod(fmmodulated_signal,fc,fs,Kfm); recoverpm = fmdemod(pmmodulated_signal,fc,fs,Kfm); %计算信噪比 signal_power_recover = mean(abs(original_signal).^2); noise_power_fm_recover = mean(abs(recoverfm - original_signal).^2); snr_fm_db = 10 * log10(signal_power_recover / noise_power_fm_recover); noise_power_pm_recover = mean(abs(recoverpm - original_signal).^2); snr_pm_db = 10 * log10(signal_power_recover / noise_power_pm_recover); %画出恢复信号时域图 figure; plot(t, recoverfm); xlabel('Time (s)'); ylabel('Amplitude'); title('recovered signal 时域图 fm'); figure; plot(t, recoverpm); xlabel('Time (s)'); ylabel('Amplitude'); title('recovered signal 时域图 pm'); figure; plot(t,[original_signal recoverfm recoverpm]); legend('Original signal','Recovered fmmod signal','Recovered pmmod signal'); xlabel('Time (s)') ylabel('Amplitude') %画出恢复信号频谱图 figure;plot(f,[20*log10(abs(fftshift(fft(original_signal)))) 20*log10(abs(fftshift(fft(recoverfm))))]); legend('Original signal','Recovered fmmod signal'); xlabel('Frequency (Hz)'); ylabel('Amplitude (dB)'); % % audiowrite('recovered_signal_fm.wav',recoverfm,fs); % audiowrite('recovered_signal_pm.wav',recoverpm,fs);