mri_model.m

This is GMS down upper converter and down converter in simulink. you may understand the structure in

% 
%  $Id: mri_model.m,v 1.1.2.3 2008/06/12 16:12:55 igork Exp $
% 
%  DISCLAIMER OF LIABILITY
%  
%  This text/file contains proprietary, confidential
%  information of Xilinx, Inc., is distributed under license
%  from Xilinx, Inc., and may be used, copied and/or
%  disclosed only pursuant to the terms of a valid license
%  agreement with Xilinx, Inc. Xilinx hereby grants you a 
%  license to use this text/file solely for design, simulation, 
%  implementation and creation of design files limited 
%  to Xilinx devices or technologies. Use with non-Xilinx 
%  devices or technologies is expressly prohibited and 
%  immediately terminates your license unless covered by
%  a separate agreement.
% 
%  Xilinx is providing this design, code, or information 
%  "as-is" solely for use in developing programs and 
%  solutions for Xilinx devices, with no obligation on the 
%  part of Xilinx to provide support. By providing this design, 
%  code, or information as one possible implementation of 
%  this feature, application or standard, Xilinx is making no 
%  representation that this implementation is free from any 
%  claims of infringement. You are responsible for 
%  obtaining any rights you may require for your implementation. 
%  Xilinx expressly disclaims any warranty whatsoever with 
%  respect to the adequacy of the implementation, including 
%  but not limited to any warranties or representations that this
%  implementation is free from claims of infringement, implied 
%  warranties of merchantability or fitness for a particular 
%  purpose.
% 
%  Xilinx products are not intended for use in life support
%  appliances, devices, or systems. Use in such applications is
%  expressly prohibited.
% 
%  Any modifications that are made to the Source Code are 
%  done at the user抯 sole risk and will be unsupported.
% 
% 
%  Copyright (c) 2008 Xilinx, Inc. All rights reserved.
% 
%  This copyright and support notice must be retained as part 
%  of this text at all times. 

%% mri_model.m
% This is a model for the multi-carrier(MC) GSM DDC (Digital DownConverter).
% The model attempts to replicate the hardware as much as possible, by 
% quantizing appropriately.

%% Set modelling options
Quantize = 1;     % double (0) or quantized (1) coefficients and data 
ShowPlots = 4;    % Show filter response and PSD plots
                  % 0 = None
                  % 1 = Essential (input, output, etc.)
                  % 2 = Important (selected filter responses)
                  % 3 = All (include filter design exploration)

use_cfir = 1;

%% Set up the design parameters
fadc =  170*1.0e6;            % Sample rate of ADC
flarmor = 67e6;               % Larmor frequency
fwidth  = 2e5;                % Spectral width of Larmor frequency
input_adc_bitwidth = 14;      
available_input_pins = 300;   % number of pins on FPGA allocated for ADC

%% Testbench signal parameters 
%  These parameters don't affect the MRI DDC structure but may
%  affect the accuracy of the measurements

freq_res = 64;                % Frequency resolution
max_sample_length = 2^18;     % Sample length, recommended maximum
signal_to_noise = 5e2;        % Input signal-to-noise ratio


sample_length = 2^floor(log2(fadc*freq_res/fwidth))   % Calculated 
if (sample_length > max_sample_length)
    sample_length = max_sample_length
end

%% Calculated model parameters

internal_bitwidth = input_adc_bitwidth+3;    % Internal bitwidth, used after the first multiplication
internal_fixed_point = input_adc_bitwidth+2; % Fixed point position
ast = 20*log10(2^internal_bitwidth); % Required stopband attenuation 

fclk = fadc;                  % Clock rate is equivalent to the ADC sampling rate
m_cfir = 2;                   % Rate change of CIC Compensation FIR

% Total rate change of DDC is calcultated with the assumption 
% that the signal should take a quarter of the first CIC response lobe
% (this might be changed if lower ripple level is required)
m_ddc = 2^floor(log2(fadc/((fwidth*m_cfir)*4)));   

m_cic = m_ddc/m_cfir;         % Rate change of CIC
fs_cic  = fadc;               % Input sampling frequency of CIC
fs_cfir = fs_cic/m_cic;       % Input sampling frequency of CFIR
      
% number of input channels is calculated taking into account the number of
% pins available and rate downsampling in CIC
channel = min(max(floor(m_cic/16),1), floor(available_input_pins/input_adc_bitwidth));
                        

%% Run filter design
[ hcic_ddc hcfir_ddc] = mri_filter_model(...
    flarmor, fwidth,...
    fclk, fs_cic, m_cic,...
    use_cfir,fs_cfir,m_cfir,...
    channel, ast, ShowPlots);


%% Generate input
clear sample;
clear harm;

sample(1:sample_length)=0;
harm(1:sample_length)=0;


% Building a house-shaped spectrum signal, with the spectrum amplitude of
% the shape:
%   ___ /\
%   | |/  \
%   |      \
%  /        \
% /          \____
% |               |
% |               |
% |               |
% |               |____
%
% The goal is to have an easily recognisable spectrum with flat plateaus in
% order to see the ripple effects.

for ii=1:freq_res
    if ((ii-1)<freq_res/8 || ((ii-1)<3*freq_res/8 && (ii-1)>=2*freq_res/8))
        input_amp(ii)=0.5+0.5*(ii-1)/(3*freq_res/8); % left slope
    end
    if((ii-1)<2*freq_res/8 && (ii-1)>=freq_res/8)    
        input_amp(ii) = 1;                           % chimney
    end
    if((ii-1)<6*freq_res/8 && (ii-1)>=3*freq_res/8)
        input_amp(ii) = 1-0.5*(ii-1-3*freq_res/8)/(3*freq_res/8); % right slope
    end
    if(ii>6*freq_res/8)
        input_amp(ii)=0.5;                           % extension
    end
    
end

% Add random phase mix
input_phase= rand(1, freq_res)*2*pi;
j=sqrt(-1); 
period = 1/fadc;
freq_step = 1/(period*sample_length);


%And do inverse FFT
for i=1 : sample_length
    freq=freq_step*(i-1);
    if(freq>=flarmor)
        if (freq<flarmor+fwidth)
            if(ceil((freq-flarmor)/fwidth*freq_res)<freq_res)
                harm(i)=(input_amp(floor((freq-flarmor)/fwidth*freq_res)+1)*exp(j*input_phase(ceil((freq-flarmor)/fwidth*freq_res))));
                max_freq=freq;
            end
        end
    else
        harm(i)=0;
    end
end
signal_max_freq=max_freq/1e6
spectrum_max_freq=freq/1e6
sample=ifft(harm);


figure; plot(real(sample));
sample=real(0.98*sample/max(abs(sample)));
figure; plot(real(sample));
test = fft(sample);
figure; plot(abs(test));

ip_quant = quantizer('fixed','round','wrap',[input_adc_bitwidth,input_adc_bitwidth-1]);
sample = quantize(ip_quant,sample);
clear mc_mixed_iq;
mc_mixed_iq_clean=[sample, sample];


step_delay=17; % Pipeline delay in hardware

clear input_signal;
for chan_index=1:channel
    if(chan_index==1)
        mc_mixed_iq_chan=mc_mixed_iq_clean+(rand(1, length(mc_mixed_iq_clean))-0.5)*max(abs(mc_mixed_iq_clean))/signal_to_noise;
        mc_mixed_iq=mc_mixed_iq_clean+(rand(1, length(mc_mixed_iq_clean))-0.5)*max(abs(mc_mixed_iq_clean))/signal_to_noise;
    else
       mc_mixed_iq_chan=(rand(1, length(mc_mixed_iq_clean))-0.5)*max(abs(mc_mixed_iq_clean))/signal_to_noise;
    end
    input_signal_single_channel=[[0:length(mc_mixed_iq_chan)-1]; real(mc_mixed_iq_chan)]';
    input_signal{chan_index}=input_signal_single_channel;
end




%% DDS
dds_quant = quantizer('fixed','floor','wrap',[internal_bitwidth,internal_bitwidth-1]);
j=sqrt(-1);
ph = ([0:length(mc_mixed_iq)-1])*floor((flarmor)/freq_step)*freq_step/fadc;
sinusoid = real((exp(-j.*2.*pi.*ph).*1.0))/2; % Cosine only, divided by 2 to match DDS
%sinusoid = (exp(-j.*2.*pi.*ph).*1.0);        % Complex multiplication
sinusoid = quantize(dds_quant, sinusoid);
sin_calc=[[0:length(mc_mixed_iq)-1]; sinusoid]';
if (exist('dds_cosine'))
  sinusoid=dds_cosine(step_delay:length(sinusoid)+step_delay-1)'; %get the sinusoid from the model
end

%% Complex Mixer

% Mutiply the iq data matrix by the sinusoid matrix
mixed_iq = mc_mixed_iq.*sinusoid;

mult_quant = quantizer('fixed','floor','wrap',[internal_bitwidth, internal_bitwidth-1]);

mixed_iq = quantize(mult_quant,mixed_iq);

% Debug
% if (exist('mult_out_hw'))
%     printf('Difference after the mixer: %f\n' ,max(abs(mult_out_hw-mixed_iq))); 
% end
%     
    

%% CIC
cic_quant=quantizer('fixed','floor','wrap',[internal_bitwidth,internal_fixed_point]);
ycic = filter(hcic_ddc,mixed_iq);
ycic_norm = double(ycic)/gain(hcic_ddc);
ycic_norm=quantize(cic_quant,ycic_norm);

% Visualise result
ycic_length=length(ycic_norm);
fft_ycic_norm = fft(double(ycic_norm(ycic_length/4+1:3*ycic_length/4)));
figure; plot(abs(fft_ycic_norm)); title('CIC (Absolute)');
figure; plot(20*log10(abs(fft_ycic_norm)/max(abs(fft_ycic_norm)))); title('CIC (dB)');

%% CFIR
cfir_quant=quantizer('fixed','floor','wrap',[internal_bitwidth,internal_fixed_point]);
yc = filter(hcfir_ddc,ycic_norm);  
yc_quant=quantize(cfir_quant,double(yc));

% Visualise result
yc_length=length(yc);
fft_yc = fft(double(yc_quant(yc_length/4+1:3*yc_length/4)));
figure; plot(abs(fft_yc)); title('CFIR (Abs)');
figure; plot(20*log10(abs(fft_yc)/max(abs(fft_yc)))); title('CFIR (dB)');


%% Save output to file for analysis

% Save it to a file for use with the DDC model
save('mri_model_output.mat','yc_quant');