%% PSD Comparison of ADMM1, ADMM2, CFR and seperated OICF
%%% 经过SSPA后的性能比较
clc;clear;
%% System Parameter
% load system_parameter.mat K1 K2 deltaK2 B N1 N2 G L Finv1 Finv2 F1 F2 

global X1 X2 x_mixed
% load Index4CCDF_Comparison.mat Index X1 X2 x_mixed

M = 4;
modu = [0 0 1 1; 0 1 0 1];
% load CCDF_Comparison_ADMM_P5_5000.mat zmixed_m_A1 zmixed_m_A2 
% load CCDF_Comparison_OICF_P5_5000.mat zmixed_m_O
% load CCDF_Comparison_ICF_P5_5000.mat zmixed_m_R
% save PSD_comparison.mat...
%     x_mixed zmixed_m_A1 zmixed_m_A2 zmixed_m_R zmixed_m_O
% load PSD_comparison.mat...
%     x_mixed zmixed_m_A1 zmixed_m_A2 zmixed_m_R zmixed_m_O
load(['TXPackets\rx_for_User' num2str(1) '.mat'],'rx'); % 加载用户数据
x_mixed = rx;
symNum = length(x_mixed);
%% 取出PAPR最大的100个数据
[PAPR_dex,MeanPower_Orignal] = calculate_PAPR(x_mixed);
symNum = 1000;
for i = 1:symNum
    [a(i),b(i)] = max(PAPR_dex);
    PAPR_dex(b(i)) = 0;
end
% b = randi([1 5000],1,symNum);
x_mixed = x_mixed(:,b);[PAPR_dex,~] = calculate_PAPR(x_mixed);10*log10(mean(PAPR_dex));
% zmixed_m_A1=zmixed_m_A1(:,b);[PAPR_dex,~] = calculate_PAPR(zmixed_m_A1);10*log10(mean(PAPR_dex))
% zmixed_m_A2=zmixed_m_A2(:,b);[PAPR_dex,~] = calculate_PAPR(zmixed_m_A2);10*log10(mean(PAPR_dex))
% zmixed_m_R=zmixed_m_R(:,b); [PAPR_dex,~] = calculate_PAPR(zmixed_m_R);10*log10(mean(PAPR_dex))
% zmixed_m_O=zmixed_m_O(:,b); [PAPR_dex,~] = calculate_PAPR(zmixed_m_O);10*log10(mean(PAPR_dex))
%% SSPA Model
p = 3;IBO = 5;L=4;
% Original symbol with SSPA
x_origin_out = SSPA_model(p,IBO,x_mixed);
[~,MeanPower] = calculate_PAPR(x_origin_out);
x_origin_out = x_origin_out./sqrt(MeanPower*L);
% % ADMM1 symbol with SSPA
% x_A1_out = SSPA_model(p,IBO,zmixed_m_A1);
% [~,MeanPower] = calculate_PAPR(x_A1_out);
% x_A1_out = x_A1_out./sqrt(MeanPower*L);
% % ADMM2 symbol with SSPA
% x_A2_out = SSPA_model(p,IBO,zmixed_m_A2);
% [~,MeanPower] = calculate_PAPR(x_A2_out);
% x_A2_out = x_A2_out./sqrt(MeanPower*L);
% % CFR symbol with SSPA
% x_R_out = SSPA_model(p,IBO,zmixed_m_R);
% [~,MeanPower] = calculate_PAPR(x_R_out);
% x_R_out = x_R_out./sqrt(MeanPower*L);
% % OICF symbol with SSPA
% x_O_out = SSPA_model(p,IBO,zmixed_m_O);
% [~,MeanPower] = calculate_PAPR(x_O_out);
% x_O_out = x_O_out./sqrt(MeanPower*L);
%% PSD plot
mode = 1; % 1.周期图法：periodogram 2.周期图法加汉明窗 2.文氏窗法：pwelch 
N = (512+36)*2;shift = -N/8;

[pxx,freq] = sprectrum(x_origin_out,mode,N);
pxx = circshift(pxx,shift);
plot(freq/pi,10*log10(pxx),'k','LineWidth',1);hold on;grid on

% [pxx,freq] = sprectrum(x_A1_out,mode,N);
% pxx = circshift(pxx,shift);
% plot(freq/pi,10*log10(pxx),'r','LineWidth',1);
% 
% [pxx,freq] = sprectrum(x_A2_out,mode,N);
% pxx = circshift(pxx,shift);
% plot(freq/pi,10*log10(pxx),'b','LineWidth',1);
% 
% [pxx,freq] = sprectrum(x_R_out,mode,N);
% pxx = circshift(pxx,shift);
% plot(freq/pi,10*log10(pxx),'.-r','LineWidth',1);
% 
% [pxx,freq] = sprectrum(x_O_out,mode,N);
% pxx = circshift(pxx,shift);
% plot(freq/pi,10*log10(pxx),'g','LineWidth',1);
% 
% xlim([-1 1]);ylim([-38 0.5])
% legend('Origin','ADMM1','ADMM2','CFR','OICF')
% xlabel('Normalized Frequency')
% ylabel('PSD (dB)')

%% APPENDIX FUNCTIONS
function [x_out] = SSPA_model(p,IBO,x)
    [~,Mean_Power] = calculate_PAPR(x); 
%     PAPR_out_O = 10*log10(PAPR_dex);
    C = sqrt(Mean_Power)*10^(IBO/20);
    x_out = x./ (1+(abs(x)./C).^(2*p)).^(1/(2*p));
end

function [pxx,freq] = sprectrum(x,mode,N)
    switch mode 
        case 1
            [pxx,freq] = periodogram(x *sqrt(2*pi),[],N,'centered');
        case 2
            [pxx,freq] = periodogram(x *sqrt(2*pi),hamming(length(x(:,1))),N,'centered');
        case 3
            [pxx,freq] = pwelch(x *sqrt(2*pi),[],[],N,'centered');
    end
    pxx = mean(pxx,2);
end

function [PAPR_dex,Mean_Power] = calculate_PAPR(x)
% PAPR_dex 计算得到所有符号各自的PAPR
% Mean_Power 符号的平均功率
    Signal_Power = abs(x).^2;
    Peak_Power   = max(Signal_Power,[],1); % norm(x,Inf)^2
    Mean_Power   = mean(Signal_Power,1);   % norm(x,2)^2/(K*IF)    
    PAPR_dex     = Peak_Power./Mean_Power;
end

function [EVM1,EVM2,EVM3,EVM_sum] = calculate_EVM(s1,s2,s3,symNum)
    % 计算得到所有符号各自的EVM
    global X1 X2_1 X2_2
    EVM1=zeros(1,symNum);EVM2=zeros(1,symNum);EVM3=zeros(1,symNum);EVM_sum=zeros(1,symNum);
    for i = 1:symNum
        EVM1(i) = norm(X1(:,i)-s1(:,i),2)/norm(X1(:,i),2);
        EVM2(i) = norm(X2_1(:,i)-s2(:,i),2)/norm(X2_1(:,i),2);
        EVM3(i) = norm(X2_2(:,i)-s3(:,i),2)/norm(X2_2(:,i),2);
        EVM_sum(i) = EVM1(i) + EVM2(i) + EVM3(i);
    end
end