title: Problem in undersampling M = 4 in Hybrid reconstruction in MRI parallel imaging in MATLAB
body: <p>I am learning image reconstruction in MRI parallel imaging. I am doing programing in Maltab about hybrid reconstruction on the paper &quot;A generalized approach to parallel magnetic resonance imaging (D.K. Sodickson, 2001)&quot;.
  My code only works on undersampling M=2, but it does not work on M=4. I do not know why.
  Could you help me?</p>
  <pre><code>clear;clc
  Nc=8;FOV=256;Nx=FOV;Ny=FOV;
  ph=phantom('modified shepp-logan',FOV);
  c_sens=zeros(Ny,Nx,Nc);
  for n=1:Nc
  for i=1:Ny
      for j=1:Nx
          mean=(FOV/(Nc+1))*n;
          var=30;
          a = 1/(var*(2*pi)^(0.5));
          b = (i-mean)^2;
          d = 2*((var)^2);
          k = (-1)*(b/d);
          c_sens(i,j,n)=a*exp(k);
      end
  end
  end
  c_img=zeros(Ny,Nx,Nc);
  for n=1:Nc
  c_img(:,:,n) = ph.*c_sens(:,:,n); 
  end
  raw=fftshift(fft2(fftshift(c_img)));
  %-----------------------------------------------------------------
  M=4;
  sampling_location = (1:M:Ny); % Ny/2+1 is sampled.
  reduced_k_space = zeros(Ny/M, Nx, Nc);
  for n = 1:Nc
  reduced_k_space(:, :, n)=raw(sampling_location,:, n);
  end
  block_size=2;
  total_block = Nc*block_size;        %total number of subblocks
  exps = zeros(block_size, Ny);
  y = -128:127;
  for m = 1:block_size
  exps(m,:) = exp(-2*pi*1i*(m-1)*y/Ny);
  end
  S_comp = zeros(Nx,Ny);
  for x=1:Nx
  C = [];
  for n=1:Nc
      C = [C; c_sens(:,x,n)'];    
  end
      
  B = [];
  tmp_b = [];
  
  %for ky = 0:M:Ny
  for ky = 0:M:(Ny/M) - 1
      for i=1:Nc
          temp_B = C(i,:) .* exp(-2*pi*1i*ky*y/Ny); % coil sensitivity * gradient modulation 
          tmp_b = [tmp_b; temp_B];
      end
      B = [B; tmp_b];
      tmp_b=[];
  end
  
  for h=1:total_block/2
      b_comb(h,:) = B(h,:);
  end 
  
  % Weights
  weights = exps*pinv(b_comb);
  
  permuted_raw = permute(reduced_k_space, [3, 1, 2]);
  reshaped_raw = reshape(permuted_raw, Nc*Nx/M, Ny);
          
  for n=1:(Ny/block_size)
      S_comp(block_size*(n-1)+1 : block_size*n, x) = weights * reshaped_raw((total_block/2)*(n-1) + 1 : (total_block/2)*n, x);
  end
  end
  
  % final reconstructed image
  reconimg = ifftshift(ifft2(ifftshift(S_comp)));
  figure(1)
  imshow(abs(reconimg),[])
  </code></pre>