MIMO-OFDM系统中随着SNR增加系统的频谱效率仿真,包含DFT码本,beam训练以及波束扫描-源码

clc; clear; close all; warning off; addpath(genpath(pwd)); %% System parameters fc = 30e9; % carrier frequency N = 64; % number of subcarriers L = 20; % number of clusters M = 20; % number of non-resolvable paths per cluster numRF = 1; % 1 data stream % Antenna arrays - use isotropic antenna elements % Transmitter Ntv = 8; Nth = 8; Nt = Ntv*Nth; % Define tapers to reduce sidelobes dBdown = 30; % dB taperz = chebwin(Ntv,dBdown); tapery = chebwin(Nth,dBdown); tap = taperz*tapery.'; % Multiply vector tapers to get 8-by-8 taper values arrayTx = phased.URA('Size',[Ntv Nth],'ElementSpacing',[0.5*physconst('LightSpeed')/fc 0.5*physconst('LightSpeed')/fc],'Taper',tap); posTx = getElementPosition(arrayTx); % Receiver Nrv = 4; Nrh = 4; Nr = Nrv*Nrh; % Define tapers to reduce sidelobes dBdown = 30; % dB taperz = chebwin(Nrv,dBdown); tapery = chebwin(Nrh,dBdown); tap = taperz*tapery.'; % Multiply vector tapers to get 8-by-8 taper values arrayRx = phased.URA('Size',[Nrv Nrh],'ElementSpacing',[0.5*physconst('LightSpeed')/fc 0.5*physconst('LightSpeed')/fc],'Taper',tap); posRx = getElementPosition(arrayRx); %% DFT codebook eleSpacing = 0.5; % element spacing, normalized by wavelength [beamTx,beamAngleTx,beamAngleElTx,beamAngleAzTx,beamElTx,beamAzTx] = getDFTCodebook(Ntv,Nth,eleSpacing,eleSpacing); [beamRx,beamAngleRx,beamAngleElRx,beamAngleAzRx,beamElRx,beamAzRx] = getDFTCodebook(Nrv,Nrh,eleSpacing,eleSpacing); %% Perform beam training % Generate MIMO channels [h,arrayResponseTx,arrayResponseRx,pathGain] = MIMOChan(Nt,Nr,L,M,posTx,posRx); H = 1/sqrt(N)*fft(h,N,3); % frequency-domain channel % Beam sweeping rPower = zeros(Nt,Nr,N); for sc = 1:N tempH = H(:,:,sc); % Nr x Nt for tb = 1:Nt % search all beam pairs f = beamTx(:,tb); % Nt x 1 for rb = 1:Nr w = beamRx(:,rb); % Nr x 1 rPower(tb,rb,sc) = sum(abs(w'*tempH*f).^2); % sum over all RF chains end end end rPowerAve = mean(rPower,3); % average over subcarriers power = rPowerAve; % Order beam pairs in descending order of receive power beamTable = zeros(Nr*Nt,2); tbIdxVec = 1:Nt; % transmit beam index rbIdxVec = 1:Nr; % receive beam index for bp = 1:(Nt*Nr) [tB,rB] = find(power == max(max(power))); power(tB(1),rB(1)) = -Inf; beamTable(bp,1) = tbIdxVec(tB(1)); beamTable(bp,2) = rbIdxVec(rB(1)); end beamSelected = beamTable(1,:); %% Data rate SNR = 0:2:10; SNR_linear = 10.^(SNR./10); dataRate = zeros(N,numel(SNR)); f = beamTx(:,beamSelected(:,1).'); w = beamRx(:,beamSelected(:,2).'); for snr = 1:numel(SNR) for sc = 1:N tempH = H(:,:,sc); dataRate(sc,snr) = abs(log2(det(eye(numRF)+SNR_linear(snr)/numRF*w'*tempH*f*f'*tempH'*w))); end end dataRate = mean(dataRate,1); figure(); plot(SNR,dataRate,'b-o'); grid on; xlabel('SNR (dB)'); ylabel('Spectral efficiency (bit/s/Hz)');