16APSK_16APSK仿真文件_16apsk_

clear;close all;clc; fs = 4.096e6; %采样率 T = 1/fs; %采样间隔 M = 16; SNR = 100; %信噪比 Rbit = 16e3; %比特速率 Rs = Rbit/log2(M); %码元（符号）速率 Ts = 1/Rs; dataNum = 320; len = dataNum*fs*Ts; N_samp = fs*Ts; t = (0:len - 1)*T;%时间 A = 20; %消息幅度 Sm = (randi(M,dataNum,1)-1)'; %消息 % Sm = [0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15]; MM = [4 12]; RAD = [1 2]; % apsk = apskmod(Sm,MM,RAD); constellation = [0.707106781186548 + 0.707106781186548i,-0.707106781186548 + 0.707106781186548i,0.707106781186547 - 0.707106781186548i,-0.707106781186548 - 0.707106781186548i,1.93185165257814 + 0.517638090205042i,-1.93185165257814 + 0.517638090205041i,1.93185165257814 - 0.517638090205043i,-1.93185165257814 - 0.517638090205042i,0.517638090205042 + 1.93185165257814i,-0.517638090205041 + 1.93185165257814i,0.517638090205041 - 1.93185165257814i,-0.517638090205043 - 1.93185165257814i,1.41421356237310 + 1.41421356237310i,-1.41421356237310 + 1.41421356237310i,1.41421356237309 - 1.41421356237310i,-1.41421356237310 - 1.41421356237310i]; SmMap = constellation(Sm+1); S_baseband = zeros(1,length(Sm)*N_samp); for i = 1:length(Sm) S_baseband(N_samp*(i-1)+1) = SmMap(i); end s_out = apskdemod(SmMap,MM,RAD); fc = 40000; %载波频率 %figure;plot(exp(1i*2*pi*fc*t));hold on;plot(exp(1i*2*pi*fc*t+pi/2)); c_base = exp(1i*2*pi*fc*t); phi = zeros(1,length(Sm)); for m = 1:length(Sm) phi(m) = 2*pi*Sm(m)/M; c(N_samp*(m-1)+1:N_samp*m) = exp(1i*phi(m)); end psk = A*c_base.*c; %psk信号 % psk = pskmod(Sm_expand,M,phi(1)); % s_out = pskdemod(psk,M,phi(1)); N_fft = 2^floor(log2(len)); psk_fft = fft(psk,N_fft); f = -fs/2:fs/N_fft:fs/2-fs/N_fft; %频谱横坐标 figure('NumberTitle','off','Name',['消息与',num2str(M),'PSK']); subplot(311) plot(Sm_expand); xlabel('数据长度') ylabel('幅度') set(gca,'FontSize',18) subplot(312) plot(real(psk),'.-b'); hold on; plot(imag(psk),'.-r'); plot(sqrt(real(psk).^2 + imag(psk).^2),'.-g'); legend('实部','虚部','包络') xlabel('数据长度') ylabel('幅度') title([num2str(M),'PSK信号实部虚部及包络']); set(gca,'FontSize',18) subplot(313) figure('NumberTitle','off','Name','消息频谱'); plot(f,fftshift(abs(psk_fft)).^4); xlabel('频率/Hz') ylabel('幅度') % title([num2str(M),'PSK信号频谱']); set(gca,'FontSize',18) %% 功率计算与加噪 Pow_S = sum(abs(psk).^2)/length(psk); %信号平均功率 Pow_noise = 1/10^(SNR/10)*Pow_S; %噪声功率 S_noise = sqrt(Pow_noise*0.5)*randn(1,length(psk)) + 1i*sqrt(Pow_noise*0.5)*randn(1,length(psk)); S_Rx = psk + S_noise; N_fft = 2^floor(log2(len)); psk_fft = fft(S_Rx,N_fft); f = -fs/2:fs/N_fft:fs/2-fs/N_fft; %频谱横坐标 figure('NumberTitle','off','Name','消息频谱'); plot(f,fftshift(abs(psk_fft))); xlabel('频率/Hz') ylabel('幅度') %% 变频 c_opposite = exp(1i*2*pi*(-fc+5e6)*t); S_convf = S_Rx.*c_opposite; %% 滤波 Rp_LPF = 1; Rs_LPF = 80; fpass_LPF = 0.4e6; %1.2e6 fcut_LPF = 0.5e6; %1.5e6 [n,fo,mo,w] = firpmord([fpass_LPF,fcut_LPF],[1,0],[(10^(Rp_LPF/20)-1)/(10^(Rp_LPF/20)+1),10^(-Rs_LPF/20)],fs); b = firpm(n,fo,mo,w); % keyboard yy = zeros(1,length(b)); S_convf_zeros = [S_convf,yy]; S_filt_zeros = filter(b,1,S_convf_zeros); if(mod(length(b),2)==0) S_filt = S_filt_zeros(:,length(b)/2+1:length(S_convf)+length(b)/2); else S_filt = S_filt_zeros(:,(length(b)+1)/2+1:length(S_convf)+(length(b)+1)/2); end %% 解调 ang_IQ = angle(S_filt); out = ang_IQ.*M/2/pi; oliugei = round(out); oliugei_zheng = zeros(1,length(oliugei)); for o = 1:length(oliugei) if oliugei(o)<0 oliugei_zheng(o) = oliugei(o) + M; else oliugei_zheng(o) = oliugei(o); end end S_filt_averary = zeros(1,length(Sm)); for in = 1:length(Sm) S_filt_averary(in) = 1/N_samp*sum(oliugei_zheng(N_samp*(in-1)+1:N_samp*in)); end %取包络 %S_abs = abs(S_filt_averary); %四舍五入取整 s_out = round(S_filt_averary); %--------- %% 变频前后星座图 figure('NumberTitle','off','Name','变频前后星座图'); subplot(211) plot(S_convf/max(abs(S_convf)),'*'); hold on; % plot(c,'r'); title('信号星座图'); legend('变频后信号','原始消息') set(gca,'FontSize',18) subplot(212) plot(S_Rx/max(abs(S_Rx)),'.'); title([num2str(M),'PSK变频前归一化星座图']); set(gca,'FontSize',18) %% 滤波 % Rp_LPF = 1; % Rs_LPF = 80; % fpass_LPF = 1.3e6; % fcut_LPF = 1.5e6; % [n,fo,mo,w] = firpmord([fpass_LPF,fcut_LPF],[1,0],[(10^(Rp_LPF/20)-1)/(10^(Rp_LPF/20)+1),10^(-Rs_LPF/20)],fs); % b = firpm(n,fo,mo,w); % % yy = zeros(1,length(b)); % S_convf_zeros = [S_convf,yy]; % % S_filt_zeros = filter(b,1,S_convf_zeros); % if(mod(length(b),2)==0) % S_filt = S_filt_zeros(:,length(b)/2+1:length(S_convf)+length(b)/2); % else % S_filt = S_filt_zeros(:,(length(b)+1)/2+1:length(S_convf)+(length(b)+1)/2); % end S_filt = S_convf; %不用滤波器的话就用这句 %% 瞬时相位瞬时频率 ang = angle(S_convf); ang_u = zeros(1,length(ang)); for k = 1:length(ang)-1 if ang(k+1)-ang(k)>pi ang_u(k+1) = ang_u(k)-2*pi; elseif ang(k)-ang(k+1)>pi ang_u(k+1) = ang_u(k)+2*pi; else ang_u(k+1) = ang_u(k); end end ang_sum = ang+ang_u; for i = 1:length(ang_sum)-1 omiga(i+1) = (ang_sum(i+1) - ang_sum(i))*fs; end frequency = omiga/(2*pi); figure('NumberTitle','off','Name',[num2str(M),'PSK信号和瞬时频率']); subplot(221) plot(t,real(psk),'.-b'); hold on; plot(t,imag(psk),'.-r'); plot(t,sqrt(real(psk).^2 + imag(psk).^2),'.-g'); legend('实部','虚部','包络') xlabel('时间/s') ylabel('幅度') title([num2str(M),'PSK信号实部虚部及包络']); set(gca,'FontSize',18) subplot(222) plot(t,ang_sum); xlabel('时间/s') ylabel('角度') title('瞬时相位') set(gca,'FontSize',18) subplot(223) plot(t,frequency); xlabel('时间/s') ylabel('频率/Hz') title('瞬时频率'); set(gca,'FontSize',18) subplot(224) plot(f,fftshift(abs(psk_fft))); xlabel('频率/Hz') ylabel('幅度') title([num2str(M),'PSK信号频谱']); set(gca,'FontSize',18) %% 解调前后结果及误码率 %signal [num_bit_err,bit_err] = biterr(Sm,s_out,log2(M)); [num_sym_err,sym_err] = symerr(Sm,s_out); figure('NumberTitle','off','Name','解调前后结果'); stem(Sm,'*-b'); hold on; stem(s_out,'o-r'); xlabel('数据长度') ylabel('幅度') legend('原始符号','解调符号') title([num2str(M),'PSK解调前后结果','误比特率',' = ',num2str(bit_err*100,3),'%',sprintf('\n'),'误码率',' = ',num2str(sym_err*100,3),'%']) set(gca,'FontSize',18)