MATLAB Implementation
MATLAB Script
%% MEDICAL IMAGING SYSTEM IDENTIFICATION & RECONSTRUCTION
% Project: Pilot-Based CIR Estimation for Medical Image Denoising
% Methodology: 32x32 Phantom Calibration -> 128x128 Diagnostic Restoration
clear; clc; close all;
%% 1. SYSTEM PARAMETERS (Industrial Standards)
target_snr_db = 30; % Signal-to-Noise Ratio
pilot_dim = 32; % Calibration Grid
data_dim = 128; % Diagnostic Grid
% Simulated Hardware Blur (Point Spread Function)
true_psf = fspecial('gaussian', [5 5], 1.2);
%% 2. CALIBRATION PHASE (32x32 Pilot)
% Use the Shepp-Logan Phantom (Industry standard for medical imaging)
pilot_clean = phantom(pilot_dim);
% Distort the Pilot (Convolution + System Noise)
pilot_blurred = imfilter(pilot_clean, true_psf, 'circular');
% Add Signal-Dependent Noise (Simulating Quantum Mottle)
noise_var = mean(pilot_blurred(:)) / (10^(target_snr_db/10));
pilot_received = pilot_blurred + sqrt(noise_var) * randn(pilot_dim);
% ESTIMATE CIR: Using Regularized Least Squares
% H_hat = (Y .* conj(X)) ./ (|X|^2 + epsilon)
P_fft = fft2(pilot_clean);
Y_fft = fft2(pilot_received);
epsilon = 0.01; % Stabilization constant
H_est = (Y_fft .* conj(P_fft)) ./ (abs(P_fft).^2 + epsilon);
% Extract Spatial CIR and apply a Window to remove noise artifacts
cir_spatial = real(ifft2(H_est));
cir_spatial = fftshift(cir_spatial);
% Apply Hamming window to "clean" the 32x32 knowledge
[W1, W2] = meshgrid(hamming(pilot_dim));
cir_knowledge = cir_spatial .* (W1 .* W2);
%% 3. DIAGNOSTIC PHASE (128x128 Patient Data)
patient_clean = phantom(data_dim);
% Apply same hardware distortion to patient scan
patient_blurred = imfilter(patient_clean, true_psf, 'circular');
d_noise_var = mean(patient_blurred(:)) / (10^(target_snr_db/10));
patient_noisy = patient_blurred + sqrt(d_noise_var) * randn(data_dim);
%% 4. ROBUST RECONSTRUCTION (The Wiener Filter)
% A. Embed 32x32 CIR knowledge into the 128x128 frequency grid
cir_padded = zeros(data_dim, data_dim);
start_idx = (data_dim - pilot_dim)/2 + 1;
end_idx = start_idx + pilot_dim - 1;
cir_padded(start_idx:end_idx, start_idx:end_idx) = cir_knowledge;
% B. Perform LMMSE (Linear Minimum Mean Square Error) Deconvolution
H_full = fft2(ifftshift(cir_padded));
I_noisy = fft2(patient_noisy);
% Optimal Wiener Constant K (Noise Power / Signal Power)
K = 1 / (10^(target_snr_db/10));
restored_fft = I_noisy .* (conj(H_full) ./ (abs(H_full).^2 + K));
patient_restored = real(ifft2(restored_fft));
%% 5. VISUALIZATION & ANALYTICS
figure('Name', 'System Performance Analysis', 'Color', 'w', 'Position', [100 100 1200 500]);
subplot(1,3,1); imshow(patient_noisy, []);
title(['Raw Scan (', num2str(target_snr_db), 'dB SNR)']);
xlabel('Input: Distortion + Quantum Mottle');
subplot(1,3,2); mesh(fftshift(abs(H_full)));
title('Learned System Response');
zlabel('Magnitude'); colormap jet;
subplot(1,3,3); imshow(patient_restored, [0 1]);
title('Reconstructed Diagnostic Image');
xlabel('Output: Sharp Decoded Result');
% Performance Metrics
psnr_val = psnr(patient_restored, patient_clean);
fprintf('--- Reconstruction Report ---\n');
fprintf('System Calibration: Successful (32x32)\n');
fprintf('Restoration Fidelity (PSNR): %.2f dB\n', psnr_val);
fprintf('Noise Suppression Ratio: %.2f%%\n', (1 - K)*100);
Output
