%SSB信号调制解调clear;clc;f0 = 1; %信源信号频率(Hz)E0 = 1; %信源信号振幅(V)E = 1; %载波分量振幅(V)fc = 10; %载波分量频率（Hz）t0 = 1; %信号时长snr = 15; %解调器输入信噪比dBdt = 0.003; %系统时域采样间隔fs = 1/dt; %系统采样频率df = 0.001; %所需的频率分辨率t = 0:dt:t0;Lt = length(t); %仿真过程中，信号长度snr_lin = 10^(snr/10);%解调器输入信噪比%-------------画出调制信号波形及频谱%产生模拟调制信号m = E*cos(2*pi*f0*t);L = min(abs(m));%包络最低点R = max(abs(m));%包络最高点%画出调制信号波形和频谱clf;figure(1);%%%画出调制信号波形subplot(411);plot(t,m(1:length(t)));axis([0,t0,-R-0.3,R+0.3]);%设置坐标范围xlabel('t');title('调制信号');set(gca,'YTick',-R:1:R);subplot(412);[M,m,df1,f] = T2F_new(m,dt,df,fs); %求出调制信号频谱[Bw_eq] = signalband(M,df,t0); %求出信号等效带宽f_start_low = fc - Bw_eq; %求出产生下边带信号的带通滤波器的起始频率f_cutoff_low = fc; %求出产生下边带信号的带通滤波器的截止频率f_start_high = fc; %求出产生上边带信号的带通滤波器的起始频率f_cutoff_high = fc + Bw_eq; %求出产生上边带信号的带通滤波器的截止频率plot(f,fftshift(abs(M))); %画出调制信号频谱%M：傅里叶变换后的频谱序列xlabel('f');title('调制信号频谱');axis([-fc-5*f0,fc+5*f0,0,max(M)+0.3]);set(gca,'XTick', -10:10:10);set(gca,'XGrid','on');%%%载波及其频谱subplot(413);c = cos(2*pi*fc*t); %载波plot(t,c);axis([0,t0,-E-0.2,E+0.2]);xlabel('t');title('载波');subplot(414); %载波频谱[C,c,df1,f] = T2F_new(c,dt,df,fs);plot(f,fftshift(abs(C))); %画出载波频谱xlabel('f');title('载波频谱');axis([-fc-5*f0,fc+5*f0,0,max(C)+0.3]);set(gca,'XTick', -10:10:10);set(gca,'XGrid','on');%%%已调信号及其频谱figure(2);subplot(321); %画已调信号u = m(1:Lt).*c(1:Lt);plot(t,u);axis([0,t0,-max(u)-0.5,max(u)+0.5]);xlabel('t');title('DSB信号');set(gca,'YTick', -max(u):1:max(u));subplot(322);[U,u,df1,f] = T2F_new(u,dt,df,fs);plot(f,fftshift(abs(U))); %画出已调信号频谱xlabel('f');title('DSB信号频谱');axis([-fc-5*f0,fc+5*f0,0,max(U)+0.3]);set(gca,'XTick', -10:10:10);set(gca,'XGrid','on');%%%滤波法产生SSB信号[H_low,f_low] = bp_f(length(u),f_start_low,f_cutoff_low,df1,fs,1);%求滤波法产生下边带需要的带通滤波器[H_high,f_high] = bp_f(length(u),f_start_high,f_cutoff_high,df1,fs,1);%[H,f] = bp_f(length(sam),f_start,f_cutoff,df1,fs,1);subplot(323);plot(f_low,fftshift(abs(H_low))); %画出带通滤波器xlabel('f');title('下边带带通滤波器');axis([-f_cutoff_low-1,f_cutoff_low+1,-0.05,1.05]);subplot(324);plot(f_high,fftshift(abs(H_high))); %画出带通滤波器xlabel('f');title('上边带带通滤波器');axis([-f_cutoff_high-1,f_cutoff_high+1,-0.05,1.05]);subplot(325);plot(f_low,fftshift(abs(H_low)));hold on; %画出带通滤波器plot(f,fftshift(abs(U))); %画出已调信号频谱axis([-fc-5,fc+5,-0.05,1.05]);xlabel('f');title('下边带信号');subplot(326);plot(f_high,fftshift(abs(H_high)));hold on; %画出带通滤波器plot(f,fftshift(abs(U))); %画出已调信号频谱axis([-fc-5,fc+5,-0.05,1.05]);xlabel('f');title('上边带信号');%%%----------------经过带通滤波器，产生单边带信号（以上边带信号为例）samuf = H_high.*U; %滤波器输出信号的频谱[samu] = F2T_new(samuf,fs); %滤波器输出信号的波形figure(3);subplot(321);plot(t,samu(1:Lt));axis([0,t0,-max(samu)-0.3,max(samu)+0.3]);xlabel('t');title('上边带信号');%%[samuf,samu,df1,f] = T2F_new(samu(1:Lt),dt,df,fs);%上边带信号频谱subplot(322);plot(f,fftshift(abs(samuf))); %画出经过理想带通滤波器后信号频谱xlabel('f');title('上边带信号频谱');axis([-fc-5*f0,fc+5*f0,0,max(samuf)+0.1]);set(gca,'XTick', -10:10:10);set(gca,'XGrid','on');%%%将已调信号送入信道%先根据所给信噪比产生高斯白噪声signal_power = power_x(samu(1:Lt)); %已调信号的平均功率noise_power = (signal_power * fs)/(snr_lin*4*Bw_eq); %求出噪声方差（噪声均值为0）noise_std = sqrt(noise_power); %噪声标准差noise = noise_std * randn(1,Lt); %产生噪声%画出信道高斯白噪声波形及频谱，此时，噪声已实现，为确知信号，可求其频谱subplot(323);plot(t,noise);axis([0,t0,-max(noise),max(noise)]);xlabel('t');title('噪声信号');subplot(324);[noisef,noise,df1,f] = T2F_new(noise,dt,df,fs); %噪声频谱plot(f,fftshift(abs(noisef))); %画出噪声频谱xlabel('f');title('噪声频谱');%%%信道中的信号%叠加了噪声的已调信号频谱sam = samu(1:Lt) + noise(1:Lt);subplot(325);plot(t,sam);axis([0,t0,-max(sam),max(sam)]);xlabel('t');title('信道中的信号');subplot(326);[samf,sam,df1,f] = T2F_new(sam,dt,df,fs); %求出叠加了噪声的已调信号频谱plot(f,fftshift(abs(samf))); %画出叠加了噪声的已调信号频谱xlabel('f');title('信道中信号的频谱');axis([-fc-5*f0,fc+5*f0,0,max(samf)+0.1]);set(gca,'XTick', -10:10:10);set(gca,'XGrid','on');%%%----------------经过带通滤波器%经过理想滤波器后的信号及其频谱DEM = H_high.*samuf; %滤波器输出信号的频谱[dem] = F2T_new(DEM,fs);%滤波器输出信号的波形figure(4);subplot(321); %经过理想带通滤波器后的信号波形plot(t,dem(1:Lt)); %画出经过理想带通滤波器后的信号波形axis([0,t0,-max(dem)-0.3,max(dem)+0.3]);xlabel('t');title('理想BPF输出信号');%%[demf,dem,df1,f] = T2F_new(dem(1:Lt),dt,df,fs);%求经过理想带通滤波器后的信号频谱subplot(322);plot(f,fftshift(abs(demf))); %画出经过理想带通滤波器后信号频谱xlabel('f');title('理想BPF输出信号频谱');axis([-fc-5*f0,fc+5*f0,0,max(demf)+0.1]);set(gca,'XTick', [-10:10:10]);set(gca,'XGrid','on');%%%--------------和本地载波相乘，即混频subplot(323);plot(t,c(1:Lt));axis([0,t0,-E-0.2,E+0.2]);xlabel('t');title('本地载波');subplot(324); %频谱载波[C,c,df1,f] = T2F_new(c(1:Lt),dt,df,fs);plot(f,fftshift(abs(C))); %画出载波频谱xlabel('f');title('本地载波频谱');axis([-fc-5*f0,fc+5*f0,0,max(C)+0.3]);set(gca,'XTick', [-10:10:10]);set(gca,'XGrid','on');%再画出混频后信号及其频谱der = dem(1:Lt).*c(1:Lt); %混频%%subplot(325); %画出混频后的信号plot(t,der);axis([0,t0,-R,R]);xlabel('t');title('混频后的信号');subplot(326);[derf,der,df1,f] = T2F_new(der,dt,df,fs); %求出混频后的信号频谱plot(f,fftshift(abs(derf))); %画出混频后的信号频谱xlabel('f');title('混频后的信号频谱');axis([-2*fc-5*f0,2*fc+5*f0,0,max(derf)+0.3]);set(gca,'XTick', [-10:10:10]);set(gca,'XGrid','on');%%%-----------------经过低通滤波器%画出理想低通滤波器figure(5);[LPF,f] = lp_f(length(der),Bw_eq,df1,fs,1); %求出低通滤波器subplot(411);plot(f,fftshift(abs(LPF))); %画出理想低通滤波器xlabel('f');title('理想LPF');axis([-f0-Bw_eq,f0+Bw_eq,-0.05,1.05]);%%%混频信号经过理想低通滤波器后的频谱及波形DM = LPF.*derf; %理想低通滤波器输出的频谱[dm] = F2T_new(DM,fs); %滤波器的输出波形subplot(412);plot(t,dm(1:Lt)); %画出经过低通滤波器后的解调波形axis([0,t0,-max(dm)-0.2,max(dm)+0.2]);xlabel('t');title('恢复信号');set(gca,'YTick', [-1:0.5:1]);set(gca,'YGrid','on');subplot(413);[dmf,dm,df1,f] = T2F_new(dm(1:Lt),dt,df,fs); %求LPF输出信号的频谱plot(f,fftshift(abs(dmf))); %画出LPF输出信号的频谱xlabel('f');title('恢复信号频谱');axis([-fc,fc,0,max(abs(dmf))+0.1]);set(gca,'XTick', [-10:10:10]);set(gca,'XGrid','on');%%subplot(414);plot(t,m(1:Lt)); %画出调制信号波形xlabel('t');title('调制信号');set(gca,'YTick', [-R:1:R]);axis([0,t0,-R-0.3,R+0.3])xlabel('t');title('调制信号');子函数%序列的傅里叶变换%各参数含义与子函数T2F中的完全相同，完成序列的傅里叶变换function [M,m,df] = fftseq(m,ts,df)fs = 1/ts;if nargin == 2n1 = 0;elsen1 = fs / df;endn2 = length(m);n = 2^(max(nextpow2(n1),nextpow2(n2)));M = fft(m,n);m = [m,zeros(1,n-n2)];df = fs / n;end%计算信号功率function p = power_x(x)%x:输入信号%p:返回信号的x功率p = (norm(x).^2)./length(x);end%将信号从频域转换到时域function [m] = F2T(M,fs)%----------------输入参数%M：信号的频谱%fs:系统采样频率%----------------输出（返回）函数%m:傅里叶逆变换后的信号，注意其长度为2的整数次幂，利用其画波形时，要注意选取m 的一部分，选取长度和所给时间序列t的长度要一致，plot(t,m(1:length(t))),否则会出错m = real(ifft(M))*fs;end%将信号从时域转换到频域function [M,m,df1,f] = T2F(m,ts,df,fs)%----------------输入参数%m：信号%ts:系统时域采样间隔、%df:所需的采样频率%fs:系统采样频率%----------------输出（返回）函数%M：傅里叶变换后的频谱序列%m:输入信号参与傅里叶变换后对应序列，需要注意的是，该序列与输入信号m的区别，其长度是不一样的，输入的m长度不一定是2的整数次幂，而傅里叶变换要求输入信号长度为2的整数次幂，%故傅里叶变换前需要对m信号进行补零操作，其长度有所增加，估输出参数中的m为补零后的输入信号，其长度与输入参数m不一样，但与M,f的长度是一样的，%并且，其与时间序列t所对应的序列m(1:length(t))与输入参数中的m是一致的。