Invalid Reconstruction using IFFT

Hello,

I am a computer science student and new to Physics and MATLAB.

I am trying to create a 1D Version of an MRI and I am unable to get a proper reconstruction IFFT. I have pasted relevant code and plots that I get. Could anyone please guide me on this.

% ========================================================================
% 1D MRI DIGITAL TWIN 
% ========================================================================
clearvars;
close all;
clc;
%% 1. SIMULATION & PLOTTING FLAGS
% ========================================================================
displayphantom     = 1;   % 1 → quick phantom plot
showWin1_FID       = 1;   % 1 → Window-1  (Real / Imag / Amp / Phase)
showWin2_Recon     = 1;   % 1 → Window-2  (Phantom / Recon / Diff)
%% 2. USER-DEFINED SCAN PARAMETERS
% ========================================================================
% --- RF Pulse Parameters ---
Alpha              = 90;      % flip-angle (deg)
RF_Duration        = 1e-3;    % RF duration  (s)
RF_Npoints         = 101;     % RF simulation samples   (1 → “instant” pulse)
RF_type            = 'rect';  % 'rect' | 'instant'
% --- Delay Parameters ---
DELAY_Duration     = 2e-3;    % first delay duration (s)
% --- Imaging & Reconstruction Parameters ---
Reconstruction_FOV = 0.20;     % imaging FOV (m)
Reconstruction_Resolution = 0.00025; % projection resolution (m)
% --- Gradient Parameters ---
ACQ_BW             = 350;     % acquisition BW (kHz), for constant gradient
%% 3. PHANTOM SETUP
% ========================================================================
Phantom_Npoints        = 2000;
Phantom_Compartments   = 5;
Compartment_Npoints    = Phantom_Npoints / Phantom_Compartments;
Magnetization          = zeros(Phantom_Npoints,1);
T1                     = zeros(Phantom_Npoints,1);   % ms
T2                     = zeros(Phantom_Npoints,1);   % ms
for idx_compartments = 1:Phantom_Compartments
    idx = (idx_compartments-1)*Compartment_Npoints + (1:Compartment_Npoints);
    switch idx_compartments
        case 1, tissue_m=0.0;  tissue_T1=600;  tissue_T2=80;
        case 2, tissue_m=0.6;  tissue_T1=1200; tissue_T2=140;
        case 3, tissue_m=1.0;  tissue_T1=1400; tissue_T2=160;
        case 4, tissue_m=0.6;  tissue_T1=800;  tissue_T2=100;
        case 5, tissue_m=0.0;  tissue_T1=1000; tissue_T2=120;
    end
    Magnetization(idx) = tissue_m;
    T1(idx) = tissue_T1;
    T2(idx) = tissue_T2;
end
signal_phantom    = Magnetization; % This represents M0, the equilibrium magnetization
Phantom_FOV_m     = 0.200; % Define phantom physical size in meters
Phantom_x         = linspace(-Phantom_FOV_m/2, Phantom_FOV_m/2, Phantom_Npoints)';
if displayphantom
    figure; plot(1000*Phantom_x, signal_phantom,'-'); grid on;
    xlabel('Position (mm)'); ylabel('M_0 (a.u.)');
    title('Original Phantom');
end

%% 4. CALCULATED PARAMETERS & CONSTANTS
% ========================================================================
% --- Physical Constants ---
Gamma = 42.58e6; % Gyromagnetic Ratio for proton (In Hz/T)
% --- Timing, Steps, and Point Calculations ---
Reconstruction_Npoints = ceil(Reconstruction_FOV / Reconstruction_Resolution);
ACQ_Npoints            = Reconstruction_Npoints;
ACQ_BW_Hz              = 1000*ACQ_BW;
ACQ_DwellTime          = 1/ACQ_BW_Hz; % Δt = 1 / BW
if strcmp(RF_type, 'instant'), RF_Npoints = 1; RF_TimeStep_s = RF_Duration; else, RF_TimeStep_s = RF_Duration/RF_Npoints; end
Delay_TimeStep_s       = ACQ_DwellTime;
Delay1_Steps           = ceil(DELAY_Duration/Delay_TimeStep_s);
% --- Angle and RF Amplitude Calculations ---
Alpha_rad              = deg2rad(Alpha);
RF_Amplitude           = Alpha_rad/(Gamma * RF_Duration); % B1 = α / (γ * T_RF)
RF_FlipAngleStep_rad   = Gamma*RF_Amplitude*RF_TimeStep_s; % Δα = γ * B1 * Δt
% --- Relaxation Factors (E1 for T1, E2 for T2) ---
T1_sec = T1/1000;
T2_sec = T2/1000;
RF_E1Relaxation   = exp(-RF_TimeStep_s ./T1_sec);    % E1 = exp(-Δt / T1)
RF_E2Relaxation   = exp(-RF_TimeStep_s ./T2_sec);    % E2 = exp(-Δt / T2)
DELAY_E1          = exp(-Delay_TimeStep_s./T1_sec);
DELAY_E2          = exp(-Delay_TimeStep_s./T2_sec);
E1_ACQ            = exp(-ACQ_DwellTime ./T1_sec);
E2_ACQ            = exp(-ACQ_DwellTime ./T2_sec);
% --- Constant Gradient & Phase Calculations ---
const_gradient = ACQ_BW_Hz/(Gamma*Reconstruction_FOV); % G = BW / (γ * FOV)
phase_inc      = 2*pi*Gamma*const_gradient*Phantom_x*ACQ_DwellTime; % ΔΦ(x) = 2π * γ * Gx * x * Δt
phase_factor   = exp(-1i*phase_inc); % Pre-calculate the phase factor e^(-iΔΦ)
%% 5. SIMULATION INITIALIZATION
% ========================================================================
voxel_Mxy = zeros(Phantom_Npoints,1);
voxel_Mz  = signal_phantom; % Start with Mz = M0
%% 6. MAIN SIMULATION LOOP
% ========================================================================
% ----- RF block ----------------------------------------------------------
for k = 1:RF_Npoints
    if strcmp(RF_type, 'instant'), dAlpha = Alpha_rad; else, dAlpha = RF_FlipAngleStep_rad; end
    if dAlpha~=0
        c=cos(dAlpha); s=sin(dAlpha);
        Mx=real(voxel_Mxy); My=imag(voxel_Mxy);
        My_new = c.*My + s.*voxel_Mz;
        Mz_new = -s.*My + c.*voxel_Mz;
        voxel_Mxy = Mx + 1i*My_new; voxel_Mz = Mz_new;
    end
    voxel_Mxy = RF_E2Relaxation.* voxel_Mxy;
    voxel_Mz  = RF_E1Relaxation.* voxel_Mz + (1-RF_E1Relaxation).*signal_phantom;
end
% ----- Delay-1 -----------------------------------------------------------
for k = 1:Delay1_Steps
    voxel_Mxy = DELAY_E2.* voxel_Mxy;
    voxel_Mz  = DELAY_E1.* voxel_Mz + (1-DELAY_E1).*signal_phantom;
end
% ----- Acquisition -------------------------------------------------------
FID = zeros(ACQ_Npoints,1);
t_acq_s = (0:ACQ_Npoints-1)'*ACQ_DwellTime;
t_acq_ms = t_acq_s*1e3;
for k = 1:ACQ_Npoints
    voxel_Mxy = phase_factor.* (E2_ACQ.* voxel_Mxy);
    voxel_Mz  = E1_ACQ.* voxel_Mz + (1-E1_ACQ).*signal_phantom;
    FID(k)=sum(voxel_Mxy);
end
%% 7. IMAGE RECONSTRUCTION
k_space = -1i * FID;
N_acquired = length(k_space);

% Perform IFFT directly on acquired k-space data without zero padding or conjugate symmetry
ProjectionFT = ifft(k_space);

% Shift and take real part for the reconstructed spatial profile
ReconstructedProfile = real(fftshift(ProjectionFT));

% Normalize the reconstructed profile to max 1 for display
if max(abs(ReconstructedProfile)) > 0
    ReconstructedProfile = ReconstructedProfile / max(abs(ReconstructedProfile));
end

pos_recon_m = linspace(-Reconstruction_FOV/2, Reconstruction_FOV/2, N_acquired);
PhantomProfile_Interpolated = interp1(Phantom_x, signal_phantom, ...
                                      pos_recon_m, 'linear', 'extrap')';
if max(abs(PhantomProfile_Interpolated)) > 0
    PhantomProfile_Interpolated = PhantomProfile_Interpolated / max(PhantomProfile_Interpolated);
end
diff_profile  = PhantomProfile_Interpolated - ReconstructedProfile;
mean_err      = mean(diff_profile);
std_err       = std(diff_profile);

%% 8. DISPLAY WINDOWS
% ========================================================================
set(0,'defaultAxesColorOrder',[0 0 0],'defaultLineLineWidth',1.2);
if showWin1_FID
    figure('Name','Window-1  FID');
    subplot(4,1,1); plot(t_acq_ms,real(FID));   ylabel('Real'); grid on; title('FID Signal (No Echo)');
    subplot(4,1,2); plot(t_acq_ms,imag(FID));   ylabel('Imag'); grid on;
    subplot(4,1,3); plot(t_acq_ms,abs(FID));    ylabel('|S|');  grid on;
    subplot(4,1,4); plot(t_acq_ms,unwrap(angle(FID)));
    ylabel('Phase'); xlabel('Time (ms)'); grid on;
end

if showWin2_Recon
    figure('Name','Window-2  Reconstruction (No Zero-Padding)');
    pos_mm = 1000*pos_recon_m;
    subplot(3,1,1); plot(pos_mm,PhantomProfile_Interpolated); grid on;
    title('Original Phantom (Interpolated & Normalized)'); ylabel('M_0 (norm)');
    subplot(3,1,2); plot(pos_mm,ReconstructedProfile); grid on;
    title('Reconstructed Profile '); ylabel('IFFT');
    subplot(3,1,3); plot(pos_mm,diff_profile); grid on;
    title(sprintf('Difference  (μ=%.3g  σ=%.3g)',mean_err,std_err));
    xlabel('Position (mm)'); ylabel('Error');
end

[dpb]

You would need to strip down to the bare minimum working example (including data so folks can duplicate it) in order for someone to try do determine what happened. It's more than most any busy volunteer is going to want to tackle to try to reconstruct from partial code without the actual data.

[Rohan]

Alright, thank you for pointing that out.