STDR.m

%% =====================================================
%  STDR Batch Simulation: Noise, Type & Distance Sweep
% =====================================================
clc; clear; close all;
model = 'STDR_3Phase_Selective';
load_system(model);

%% 1. Generate PN excitation (TX)
[pnA_vec, pnB_vec, pnC_vec, ~, Fs, ~] = generatePN3phase(1023, 20, 500e3, 'orthogonal');
Ts = 1/Fs; t = (0:length(pnA_vec)-1).' * Ts;
pnA = [t pnA_vec]; pnB = [t pnB_vec]; pnC = [t pnC_vec];
assignin('base','pnA',pnA); assignin('base','pnB',pnB); assignin('base','pnC',pnC);
v = 2e8;

%% 2. Define Experiment Matrix
faultTypes = {'LG','LL','SC','OC'};
faultDistances = [300 750 1500 2250 2700]; 
noiseLevels = [0, 1e-6, 1e-4]; 
GammaMap = containers.Map({'LG','LL','SC','OC'}, {-0.5, -0.7, -1.0, +1.0});

%% 3. Result Storage
Results = {};
runID = 1;

%% 4. Triple Nested Batch Loop
for n = 1:numel(noiseLevels)
    currentNoise = noiseLevels(n);
    
    % Update Noise blocks (Ensure these exist in Channel_3Phase)
    set_param([model '/Channel_3Phase/Noise_A'], 'Cov', num2str(currentNoise));
    set_param([model '/Channel_3Phase/Noise_B'], 'Cov', num2str(currentNoise));
    set_param([model '/Channel_3Phase/Noise_C'], 'Cov', num2str(currentNoise));
    
    for f = 1:numel(faultTypes)
        faultType = faultTypes{f};
        Gamma = GammaMap(faultType);
        
        for d = 1:numel(faultDistances)
            faultDistance = faultDistances(d);
            Ndelay = round((2 * faultDistance / v) / Ts);
            
            fprintf('Run %03d | Noise: %.1e | %s at %d m\n', runID, currentNoise, faultType, faultDistance);
            
            set_param([model '/FaultDelayConst'], 'Value', num2str(Ndelay));
            set_param([model '/Channel_3Phase/Gamma_A'],'Gain', num2str(Gamma));
            set_param([model '/Channel_3Phase/Gamma_B'],'Gain', num2str(Gamma));
            set_param([model '/Channel_3Phase/Gamma_C'],'Gain', num2str(Gamma));
            
            out = sim(model);
            
            distA = out.dist_A.Data(end); distB = out.dist_B.Data(end); distC = out.dist_C.Data(end);
            tauA  = out.tau_A.Data(end);  tauB  = out.tau_B.Data(end);  tauC  = out.tau_C.Data(end);
            
            Results(runID,:) = {currentNoise, faultType, faultDistance, ...
                                distA, distB, distC, ...
                                tauA, tauB, tauC, ...
                                abs(distA-faultDistance), abs(distB-faultDistance), abs(distC-faultDistance)};
            runID = runID + 1;
        end
    end
end

%% 5. Convert to Table
ResultsTable = cell2table(Results, 'VariableNames', ...
    {'NoiseLevel','FaultType','TrueDistance_m', ...
     'EstA','EstB','EstC', ...
     'TauA','TauB','TauC', ...
     'ErrA','ErrB','ErrC'});
writetable(ResultsTable,'STDR_Noise_Sweep_Results.csv');

%% 6. Performance Statistics (RMSE per Phase)
fprintf('\n--- Performance Summary (RMSE) ---\n');
uniqueNoise = unique(ResultsTable.NoiseLevel);
for i = 1:length(uniqueNoise)
    level = uniqueNoise(i);
    subData = ResultsTable(ResultsTable.NoiseLevel == level, :);
    rmseA = sqrt(mean(subData.ErrA.^2));
    rmseB = sqrt(mean(subData.ErrB.^2));
    rmseC = sqrt(mean(subData.ErrC.^2));
    fprintf('Noise: %.1e | RMSE (A: %.2f, B: %.2f, C: %.2f) meters\n', level, rmseA, rmseB, rmseC);
end

%% 7. Visualization (All Phases)

% --- FIGURE 1: Error Distribution across All Phases ---
figure('Color','w','Name','Error All Phases','Position', [100 100 900 700]);
subplot(3,1,1);
gscatter(ResultsTable.TrueDistance_m, ResultsTable.ErrA, ResultsTable.NoiseLevel);
ylabel('Err A (m)'); title('Phase A Error Distribution'); grid on; legend('off');
subplot(3,1,2);
gscatter(ResultsTable.TrueDistance_m, ResultsTable.ErrB, ResultsTable.NoiseLevel);
ylabel('Err B (m)'); title('Phase B Error Distribution'); grid on; legend('off');
subplot(3,1,3);
gscatter(ResultsTable.TrueDistance_m, ResultsTable.ErrC, ResultsTable.NoiseLevel);
ylabel('Err C (m)'); xlabel('True Distance (m)'); title('Phase C Error Distribution'); grid on;
legend('Ideal','Low Noise','High Noise','Location','bestoutside');
saveas(gcf, 'STDR_Error_AllPhases.png');

% --- FIGURE 2: Linearity Check (Three-Phase Overlay) ---
figure('Color','w','Name','Linearity All Phases');
hold on;
plot(ResultsTable.TrueDistance_m, ResultsTable.EstA, 'ro', 'MarkerFaceColor', 'r', 'DisplayName', 'Phase A');
plot(ResultsTable.TrueDistance_m, ResultsTable.EstB, 'gs', 'MarkerFaceColor', 'g', 'DisplayName', 'Phase B');
plot(ResultsTable.TrueDistance_m, ResultsTable.EstC, 'bv', 'MarkerFaceColor', 'b', 'DisplayName', 'Phase C');
plot(ResultsTable.TrueDistance_m, ResultsTable.TrueDistance_m, 'k--', 'LineWidth', 2, 'DisplayName', 'Ideal Line');
xlabel('True Distance (m)'); ylabel('Estimated Distance (m)');
title('STDR Estimation Consistency Across All Phases');
legend('Location','northwest'); grid on;
saveas(gcf, 'STDR_Linearity_AllPhases.png');

disp('✔ Files saved: STDR_Noise_Sweep_Results.csv, STDR_Error_AllPhases.png, and STDR_Linearity_AllPhases.png');