QPSK调制下SEFDM通信系统收发机的matlab实现

看到网上没有开源的SEFDM系统的收发机，在这里把我的matlab代码放给大家。
有关OFDM的收发机的matlab实现可以看我之前的博客：
OFDM发射机与接收机的matlab实现
接收机检测器使用的是TSVD检测器，相关论文可以看[1]

% SEFDM_TSVD

clear;
close all;

EbN0 = [0:2:30];          %Eb/N0
N_iter = 1e5;
Nframe = 3;
Nused = 16;
Nbps = 2; M = 2^Nbps;     %调制阶数 = 2/4/6; QPSK/16-QAM/64-QAM
norms = [1 sqrt(2) 0 sqrt(10) 0 sqrt(42)];  %BPSK 4-QAM 16-QAM
Target_neb = 500;    %积累了一定的错误自动停止循环
sigPow = 0;
Nfft = 16;

alpha = 0.9;
X_shift = zeros(16,1);
x_GI = zeros(Nframe*Nused,1);
Xmod_r = zeros(Nframe*Nused,1);
rho = 1;
[phi,c_ifft_fft] = phi_c_ifft_fft(Nused,rho,alpha);
xi = 15;

for i = 0:length(EbN0)
    randn('state',0);rand('state',0);
    %     初始化错误比特数和总比特数
    Neb = 0;
    Ntb = 0;
    %     迭代设计
    for m = 1:N_iter
        X = randi([0,1],1,Nused*Nframe*Nbps);   %bit:整数向量
        X = X';
        Xmod = qammod(X,M,'gray',"InputType","bit")/norms(Nbps);

        kk1 = [1:Nused];          %有效载波序号
        kk16 = 1:Nused;
        for k = 1:Nframe
           X_shift(kk16) = Xmod(kk1);
           x = phi * X_shift/rho/Nused;
           x_GI(kk1) = x;
           
           kk1 = kk1 + Nused;
           
        end
        
%         信道
        y = x_GI;
        if i == 0     %只测量信号功率
            y1 = y(1:Nframe*Nused);
            sigPow = sigPow + y1'*y1;
            continue;
        end
        snr = EbN0(i) + 10*log10(Nbps);  %式(4.28)SNR vs.Eb/N0
        noise_mag = sqrt((10.^(-snr/10))*sigPow/2);
        y_GI = y + noise_mag * (randn(size(y))+1j*randn(size(y)));

        %RX_____________________________
        kk1 = 1:Nused;
        
        for k = 1:Nframe
            [U,S,V] = svd(c_ifft_fft);
            S_xi = S;
            for i_xi = 1:length(S)
                if i_xi >xi
                    S_xi(i_xi,i_xi) = 0; 
                else
                    S_xi(i_xi,i_xi) = inv(S_xi(i_xi,i_xi));
                end
            end
            C_xi = V * S_xi * U';
            Y(kk1) = C_xi * phi'* y_GI(kk1);
            Y_shift = Y(kk1);
            Xmod_r(kk1) = Y_shift;
            
            kk1 = kk1 + Nused;
        end
%         Xmod_r = y_GI;
        X_r = qamdemod(Xmod_r*norms(Nbps),M,'gray',"OutputType","bit");
        Neb = Neb + sum(sum( X_r ~= X )); %计算误比特数
        Ntb = Ntb + Nused*Nframe * Nbps;
        if Neb > Target_neb,break;end
    end
    if i == 0
        sigPow = sigPow/Nused/Nframe/N_iter;
    else
        Ber = Neb/Ntb;
        fprintf('EbN0=%3d[dB],BER=%4d/%8d=%11.3e\n',EbN0(i),Neb,Ntb,Ber);
        Ber_buf(i) = Ber;
        if Ber < 1e-6 , break;end
    end
end

disp('Simulation finished')
EbN0dB = [0:1:30];
M = 2^Nbps;
% ber_AWGN = ber_QAM(EbN0dB,M,'AWGN');
ber_Theory = berawgn(EbN0dB, 'psk',M, 'nondiff');
semilogy(EbN0dB,ber_Theory,'r:'),hold on
semilogy(EbN0(1:i),Ber_buf,'b--s');
grid on;
axis([EbN0(1) EbN0(end) 1e-6 1])
legend('AWGN analytic','Simulation')
xlabel('EbN0[dB]'),ylabel('BER')

function [phi,c_ifft_fft] = phi_c_ifft_fft(N,rho,alpha)
    phi_row = zeros(1,N);
    phi = zeros(rho*N,N);
    for i = 1:N
       phi_row(i) =  exp(1j*2*pi*alpha*(i-1)/(rho*N));
       
    end
    for i = 1:N*rho
       phi(i,:) = phi_row.^(i-1); 
    end
    c_ifft_fft = phi'*phi/rho/N;
end

function ber = ber_QAM(EbN0dB,M,AWGN_or_Rayleigh)
N = length(EbN0dB);
sqM = sqrt(M);
a = 2*(1-power(sqM,-1))/log2(sqM);
b = 6*log2(sqM)/(M-1);
if nargin<3
   AWGN_or_Rayleigh = 'AWGN';
   
end

ber = a * Q(sqrt(b*10.^(EbN0dB/10)));

end
function y = Q(x)
    y = erfc(x/sqrt(2))/2;
end

[1]桂韬. 高频谱效率频分复用(SEFDM)在光纤通信系统中的应用[D]. 暨南大学, 2015.