Want to check the program is technically right and work accurately?

I have two signals, T_para and T_per (5*17*20) loaded from directory, I coded this program to do the following tasks:

1. Calculate the cost function
2. Apply interpolation to increase the density of a and r
3. Minimize the cost function using optimiztion method (Grid search method)
4. Plot figures with minimum cost value

Please check if there is a mistake technically or in the structure. My program code is:

clear all; clc; close all;
cloud = 'Homo';
T_para = zeros(5, 17, 20);
T_per = zeros(5, 17, 20);

% Initialize matrices
T_para_noisy = zeros(5, 17, 20);
T_per_noisy = zeros(5, 17, 20);

% Function to add AWGN
add_noise = @(signal, snr) awgn(signal, snr, 'measured');

% Define the cost function
% cost_function = @(original, noisy) sum((original(:) - noisy(:)).^2) / numel(original);
cost_function = @(original, noisy, variance) sum(((original(:) - noisy(:))./ (var(original(:)))).^2);

% Initialize arrays to store the optimal values and global minimum cost
lookup_table = [];
lookup_table_interp = [];

% Initialize global minimum cost values
global_min_cost = Inf;
optimal_a= 0;
optimal_r = 0;

% Define ranges for a and r
a_range = 0.01:0.01:0.05;
r_range = 4:20;

% Initialize cost matrix to store costs for all combinations of a and r
cost_matrix = zeros(length(a_range), length(r_range));

for h = 1000:1000:4000
    for fov = [0.2, 0.5, 1, 2, 5, 10]

        for a_idx = 1:length(a_range)
            for r_idx = 1:length(r_range)
                a = a_range(a_idx);
                r = r_range(r_idx);

                load(['./Tabledata_', cloud, '/', num2str(h), 'm-', num2str(fov), 'mrad/TABLE.mat']);

                % Add noise to the signals
                T_para_noisy = add_noise(T_para, 30); % 20dB SNR
                T_per_noisy = add_noise(T_per, 30);

                % Calculate the variance of the original signals
                var_T_para = var(T_para, 0, 'all');
                var_T_per = var(T_per, 0, 'all');

                % Calculate the cost for each pair of original and noisy signals
                cost_para = cost_function(T_para, T_para_noisy, var_T_para);
                cost_per  = cost_function(T_per, T_per_noisy, var_T_per);

                % Sum the costs
                total_cost = cost_para + cost_per;

                % Store the cost for plotting
                cost_matrix(a_idx, r_idx) = total_cost;

                % Check if this is the new global minimum
                if total_cost < global_min_cost
                    global_min_cost= total_cost;
                    optimal_a = a;
                    optimal_r = r;
                end

            end
        end
        % Store the optimal values and global minimum cost in the lookup table
        lookup_table = [lookup_table; optimal_a, optimal_r, global_min_cost];

    end
end

% Save and display the look-up table for optimal interpolated values
save('LookupTable.mat', 'lookup_table');

% Initialize global minimum cost for interpolated values
global_min_cost_interp = Inf;
optimal_a_interp = 0;
optimal_r_interp = 0;

% Define finer gridded interpolant for cost_matrix
a_interp = 0.01:0.001:0.05;
r_interp = 4:0.5:20;

% Define the original and interpolated grids
[a_grid, r_grid] = meshgrid(a_range, r_range);
[a_interp_grid, r_interp_grid] = meshgrid(a_interp, r_interp);

% Transpose the cost_matrix to match the grid dimensions
cost_matrix = cost_matrix';

% Apply gridded interpolation using spline method
cost_matrix_interp = interp2(a_grid, r_grid, cost_matrix, a_interp_grid, r_interp_grid, 'spline');

% Loop through the interpolated grid to find the global minimum
for a_interp_idx = 1:length(a_interp)
    for r_interp_idx = 1:length(r_interp)
        interpolated_cost = cost_matrix_interp(r_interp_idx, a_interp_idx); % Note the order of indices
        % Check if this is the new global minimum
        if interpolated_cost < global_min_cost_interp
            global_min_cost_interp = interpolated_cost;
            optimal_a_interp = a_interp(a_interp_idx);
            optimal_r_interp = r_interp(r_interp_idx);
        end
    end
end

% Add to the lookup table
lookup_table_interp = [lookup_table_interp; optimal_a_interp, optimal_r_interp, global_min_cost_interp];

% Save and display the look-up table for optimal interpolated values
save('LookupTable_Interp.mat', 'lookup_table_interp');

disp('Global minimum interpolated cost:');
disp(global_min_cost_interp);


% Normalize the cost_matrix_interp values
min_cost = min(cost_matrix_interp(:));
max_cost = max(cost_matrix_interp(:));
normalized_cost_matrix_interp = (cost_matrix_interp - min_cost) / (max_cost - min_cost);

%%
% Plot the normalized interpolated cost function
figure(1);
surf(a_interp_grid, r_interp_grid, cost_matrix_interp);
shading interp;
colormap parula;
colorbar;
set(gca, 'color', 'w', 'FontSize', 12, 'LineWidth', 1, 'FontWeight', 'normal');
ylabel('R_{e} interp (\mum)');
xlabel('\alpha_{e} interp (m^{-1})');
zlabel('Normalized interpolated cost function');
title('Normalized interpolated CF over \alpha_{e} interp and R_{e} interp');
set(gca, 'box', 'on', 'LineWidth', 1, 'FontName', 'Arial', 'FontSmoothing', 'on');
set(gca, 'xlim', [0.01 0.05], 'xtick', [0.01 0.02 0.03 0.04 0.05], 'ylim', [4 20], 'ytick', [4 8 12 16 20]);
grid on;
hold on;

% Create a meshgrid for plotting
[A, R] = meshgrid(a_interp, r_interp);

% Reshape the matrices to vectors for plotting
A_vector = reshape(A, [], 1);
R_vector = reshape(R, [], 1);
Re_cost_matrix_interp = reshape(cost_matrix_interp, [], 1);

% Find the global minimum in the interpolated cost matrix
[global_min_cost_interp, min_idx] = min(Re_cost_matrix_interp);
optimal_a_interp = A_vector(min_idx);
optimal_r_interp = R_vector(min_idx);

% Mark the minimum point
hold on;
surf_handle=plot3(optimal_a_interp, optimal_r_interp, global_min_cost_interp, 'ro', 'MarkerSize', 12, 'MarkerFaceColor', 'r');
view(3); % Set the viewing angle to 3D
% Bring the red point to the front
uistack(surf_handle, 'top'); % Bring the red point to the front
hold off;

% Plot the 2D Contour Plot of the interpolated cost function
figure(2);
contourf(a_interp_grid, r_interp_grid, normalized_cost_matrix_interp, 20);
colorbar;
xlabel('\alpha_{e} interp (m^{-1})');
ylabel('R_{e} interp (\mum)');
title('2D Contour Plot of interpolated normalized CF');
set(gca,'FontSize', 12, 'LineWidth', 1, 'FontWeight', 'normal');
set(gca, 'box', 'on', 'LineWidth', 1, 'FontName', 'Arial', 'FontSmoothing', 'on');
set(gca, 'xlim', [0.01 0.05], 'xtick', [0.01 0.02 0.03 0.04 0.05], 'ylim', [4 20], 'ytick', [4 8 12 16 20]);
hold on;
plot(optimal_a_interp, optimal_r_interp, 'ro', 'MarkerSize', 12, 'MarkerFaceColor', 'r'); % Mark the minimum point from interpolated grid search
% plot(optimal_a, optimal_r, 'bo', 'MarkerSize', 10, 'LineWidth', 2); % Mark the minimum point from original grid search
hold off;

% Plot the normalized cost values
figure(3);
scatter3(R_vector, A_vector, Re_cost_matrix_interp, 'filled');
xlabel('R_{e} interp (\mum)');
ylabel('\alpha_{e} interp (m^{-1})');
zlabel('Interpolated cost function');
title('Normalized interpolated CF over \alpha_{e} interp and R_{e} interp');
set(gca, 'color', 'w', 'FontSize', 12, 'LineWidth', 1, 'FontWeight', 'normal');
set(gca, 'box', 'on', 'FontName', 'Arial', 'FontSmoothing', 'on');
set(gca, 'ylim', [0.01 0.05], 'ytick', [0.01 0.02 0.03 0.04 0.05], 'xlim', [4 20], 'xtick', [4 8 12 16 20]);
grid on;

hold on;

% Set the view angle
view(60, 40); % Adjust the view angles to see the plot clearly
% Mark the minimum point with a circle
surf_handle = plot3(optimal_r_interp, optimal_a_interp, global_min_cost_interp, 'ro', 'MarkerSize', 12, 'MarkerFaceColor', 'r');

% Bring the red point to the front
uistack(surf_handle, 'top'); % Bring the red point to the front
hold off;

Answer 1

[Umar]

Hi Wiqas,

Based on the program code you provided, it seems like you are trying to perform several tasks related to signal processing and optimization. Let's break down your program structure and analyze if there are any technical or structural mistakes:

1. *Initialization and Signal Loading*: - You correctly initialize the signals `T_para` and `T_per` with zeros. - Loading signals from directories based on specific parameters (`h` and `fov`) looks good.

2. *Noise Addition and Cost Function*: - The function to add AWGN (`add_noise`) seems appropriately defined. - The cost function calculation seems valid, considering the variance normalization.

3. *Optimization using Grid Search*: - Looping over ranges of `a` and `r` to find the optimal values based on minimum cost is a standard grid search approach. - Storing optimal values in the lookup table for later analysis is a good practice.

4. *Interpolation and Global Minimum Search*: - Performing interpolation to find finer-grained optimal values is essential for improving accuracy. - The visualization of the interpolated cost function using plots is well-implemented.

5. *Final Insights*: - Ensure that the grid search ranges for `a` and `r` are appropriate for capturing the optimal values effectively. - Consider adding comments to explain the purpose of specific sections of code for better readability and maintainability.

Overall, your program seems technically sound, covering the necessary steps for calculating costs, applying interpolation, and optimizing using a grid search method. However, always ensure to validate your results thoroughly and consider edge cases during optimization.

If you encounter any issues during execution or need further assistance in optimizing your program, feel free to ask for specific guidance or debugging support. Keep up the good work!