matlab_通过数值计算得到超奈奎斯特信号处理下的信道容量

%% calc secury capacity %% change roll off factor in each slot %% Eve is noise less %% case 2: MF, no WMF clear all;close all;clc; %% parameter setting dbeta = 0.001; betamax = 0.5; betamin = 0.2; BETA = betamin:dbeta:betamax; % roll off factor, ~U(betamin,betamax) alpha = 0.8; % FTN factor % alpha = 1/(1+beta); T = 0.5; % symbol duration when no FTN % badwidth = 1/(2*T) = 1 approximatly PvsN0_dB = 0:5:30;% assum n_E = n_B % ----------------------------------------------------- % ----------------- for C_ab --------------------------- % ----------------------------------------------------- N = 1e4; % points for calc integration, in fact 2*N-1 points Kmax = 3; % freq domain shif and add number % ----------------------------------------------------- %------------------ for C_ae -------------------------- % ----------------------------------------------------- span = 50; % rRC truncat symbol number sps = 1e3; % samples per symbol = upsample rate L = 10; % ISI length -L~L %% clalc Cab, RRC pulse, FTN signling, i.i.d source w = linspace(-pi,pi,2*N-1); % -pi<= w <= pi dw = 2*pi/(2*N-1); Cabt = zeros(1,length(BETA)-1); % for calc avg for i0 = 1:length(PvsN0_dB) disp(['process PvsN0 = ',num2str(PvsN0_dB(i0))]); PvsN0 = 10^(PvsN0_dB(i0)/10); % assum n_E = n_B EsvsN0 = PvsN0*T*alpha; for i1 = 1:length(BETA)-1 beta = BETA(i1); if 0 == mod(beta,0.05) disp(['calc Cab: beta = ',num2str(beta)]); end %------------------------gen H(w) ---------------------- % dw = 2*pi*alpha*T*df; Hw = zeros(1,2*N-1); for i2 = 1:2*N-1 Gft = 0; for k = -Kmax:Kmax f = w(i2)/(2*pi*alpha*T) + k/(alpha*T); Gft = Gft + calcRcPSD(f,T,beta); end Hw(i2) = Gft; end Hw = Hw/(alpha*T); % figure; % plot(w,Hw); Cabt(i1) = sum(0.5*log2(1+2*EsvsN0*Hw))/(2*pi*alpha*T)*dw; end Cab(i0) = sum(Cabt)/(betamax - betamin)*dbeta; % bit/s end %% calc Cae, MF, no WMF g_E = rcosdesign(betamax,span,sps); % sum(g_E.^2) = 1 centerIndex = span*sps+1; Caet = zeros(1,length(BETA)-1); for i1 = 1:length(BETA)-1 % for all SNR, Cae is the same beta = BETA(i1); if 0 == mod(beta,0.05) disp(['calc Cae: beta = ',num2str(beta)]); end g_beta = rcosdesign(beta,span,sps); Rgg = conv(g_E,g_beta);% mutral correlation sampleIndex = centerIndex + round((-L:L)*alpha*sps); f = Rgg(sampleIndex); % % >>>>>>>>>>>>>>>> for test >>>>>>>> % figure; % plot(Rgg); % hold on; % stem(sampleIndex,f,'r'); % % <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< varJ = sum(f.^2) - f(L+1)^2; % no f(0) % snr = f(L+1)^2 *2*EsvsN0 / (varJ*2*EsvsN0 + 1); snr = f(L+1)^2 / (varJ); Caet(i1) = 0.5*log2(1+snr); end Cae = sum(Caet)/(betamax - betamin)*dbeta/(alpha*T); % bit/s %% clac Cs Cs = max(Cab - Cae,0); figure; plot(PvsN0_dB,Cs,'-o'); xlabel('PvsN0 dB'); ylabel('Cs bit/s'); title('alpha = 0.8, beta = [0.2 0.5], T = 0.5');