remove/correct addpath for linux comptatibility

Author: gBringout

Reco2D_IdealFFL.m

@@ -4,15 +4,23 @@
 %% 1. Define the access path to the required other packages.
 % This are typically the sphericalHarmonics and  ScannerDesign packages.
 disp('1. Define the access path to the required other packages.')
-addpath(genpath('.'))
-addpath(genpath('..\SphericalHarmonics\'))
-addpath(genpath('..\ScannerDesign\'))
+
+addpath(genpath(fullfile('.')))
+addpath(genpath(fullfile('..','SphericalHarmonics')))
+addpath(genpath(fullfile('..','ScannerDesign')))
 
 %% 2. Loading the scanner model.
 % It is made of at least the spherical harmonics coefficient, the radius in which they are define and the nominal current used for the coil.
 disp('2. Loading the scanner''s coil codel')
 load('IdealFFL.mat');
 
+% Scanner properties
+system.c11 = 2*Quadru_0.bc(1).coefficient(2,2)/Quadru_0.rhoReference*Quadru_0.current; % "Gradient" on the line
+system.aDx = Drive_X.bc(1).coefficient(1,1)*Drive_X.current;
+system.aDy = Drive_Y.bc(2).coefficient(1,1)*Drive_Y.current;
+fprintf('c11 %g T/m; aDx %g T; aDy %g T;\n',system.c11,system.aDx,system.aDx);
+fprintf('Drive peak translation x:%g m; y:%g m\n',system.aDx/system.c11,system.aDy/system.c11);
+
 %% 3. Define all the other system parameters.
 % The spatial resolution of the phantom and the system matrix. The used frequencies, time vector, sampling frequencies, noise parameters etc. This is highly dependent on the scanner you want to model.
 disp('3. Define all the other system parameters.')
@@ -124,14 +132,14 @@
 disp('5. Define the time-varying amplitude applied on the different coils.')
 c1 = cos(2*pi*system.frequencyQuadrupol*calculation.time(:).');
 c2 = sin(2*pi*system.frequencyQuadrupol*calculation.time(:).');
-c3 = cos(0.5*2*pi*system.frequencyQuadrupol*calculation.time(:).');
-c4 = sin(0.5*2*pi*system.frequencyQuadrupol*calculation.time(:).');
+c3 = sin(0.5*2*pi*system.frequencyQuadrupol*calculation.time(:).');
+c4 = -cos(0.5*2*pi*system.frequencyQuadrupol*calculation.time(:).');
 c5 = sin(2*pi*system.frequencyDrive*calculation.time(:)');
 
 c1_dt = cos(2*pi*system.frequencyQuadrupol*(calculation.time(:).'+calculation.dt));
 c2_dt = sin(2*pi*system.frequencyQuadrupol*(calculation.time(:).'+calculation.dt));
-c3_dt = cos(0.5*2*pi*system.frequencyQuadrupol*(calculation.time(:).'+calculation.dt));
-c4_dt = sin(0.5*2*pi*system.frequencyQuadrupol*(calculation.time(:).'+calculation.dt));
+c3_dt = sin(0.5*2*pi*system.frequencyQuadrupol*(calculation.time(:).'+calculation.dt));
+c4_dt = -cos(0.5*2*pi*system.frequencyQuadrupol*(calculation.time(:).'+calculation.dt));
 c5_dt = sin(2*pi*system.frequencyDrive*(calculation.time(:).'+calculation.dt));
 
 system.coefSelection_Z = ones(1,system.numberOfTimePoints);
@@ -196,7 +204,7 @@
 % 
 % figure
 % threshold = 3*10^-3;
-% for i=1:system.nbrPointPerDrivePeriod/20:system.numberOfTimePoints
+% for i=1:25:system.numberOfTimePoints
 %     image = reshape(Babs(i,:),[system.sizeXSM,system.sizeYSM]);
 %     imagesc(system.xSM,system.ySM,image);
 %     xlabel('x axis /m')

Reco2D_IdealFFP.m

@@ -4,9 +4,9 @@
 %% 1. Define the access path to the required other packages.
 % This are typically the sphericalHarmonics and  ScannerDesign packages.
 disp('1. Define the access path to the required other packages.')
-addpath(genpath('.'))
-addpath(genpath('..\SphericalHarmonics\'))
-addpath(genpath('..\ScannerDesign\'))
+addpath(genpath(fullfile('.')))
+addpath(genpath(fullfile('..','SphericalHarmonics')))
+addpath(genpath(fullfile('..','ScannerDesign')))
 
 %% 2. Loading the scanner model.
 % It is made of at least the spherical harmonics coefficient, the radius in which they are define and the nominal current used for the coil.

SignalEncoding.m

@@ -6,7 +6,7 @@
 close all
 clear all
 
-addpath('coreFunctions\')
+addpath(genpath(fullfile('.')))
 
 
 mu0         = 4*pi*1e-7; % permeability of the air.
@@ -17,18 +17,18 @@
 Ms          = 0.6/mu0; %[A/m] Magnetisation at saturation from Magnetit. From the thesis. P54
 temperature = 310;
 %signal generation and procesing
-Fs          = 25*10^6;%2.5e7; % [Hz] sampling frequency
+Fs          = 1*10^6;%2.5e7; % [Hz] sampling frequency
 fStart      = 2; % [-] staring frequency components to be displayed
-fStop       = 500; % [-] Last frequency components to be displayed. Please respect the shanon theorem :)
+fStop       = 80; % [-] Last frequency components to be displayed. Please respect the shanon theorem :)
 dt          = 1/(50*10^6); %delta used to calculate the derivative. Used to compensate the "small" sampling frequencies
 
 % MPI signal generation
 f           = 25000; % [Hz] Frequency of the signal
-nbrPeriod   = 10; % nbr of signal period
+nbrPeriod   = 4; % nbr of signal period
 time        = 0:1/Fs:nbrPeriod/f-1/Fs; % [s] time vector. Remove the last point as we start at zero.
 N           = length(time); %Number of time points
-H0          = 31831;% [A/m] Drive field amplitude (31831 A/m to generate 40 mT)
-Hoffset     = 1000;% [A/m] Offset field
+H0          = 31831/2;% [A/m] Drive field amplitude (31831 A/m to generate 40 mT)
+Hoffset     = 1000*4;% [A/m] Offset field
 B0          = H0*mu0; % [T] - H = B/mu
 
 % MPI scanner properties
@@ -40,7 +40,7 @@
 % Noise model
 kB          = 1.380650424e-23; % [J/K] Boltzmann constant 
 deltaF      = 10*10^6; % [Hz] according to Weizenecker 2007 - A simulation study...
-Rp          = 185 *10^-3; % [Ohm]  according to Weizenecker 2007 - A simulation study...
+Rp          = 0.0001*185 *10^-3; % [Ohm]  according to Weizenecker 2007 - A simulation study...
 Scaling     = 1000; % To easily scale the maximal voltage induced by the noise
 %% Particle Magnetization in function of the applied magnetic field
 
@@ -105,8 +105,8 @@
 %% Spectrum
 
 subplot(2,3,6)
-periodogram = fft(u);
-uhat = abs(2*periodogram(1:(end-1)/2+1));
+periodogram = fft(u)/size(u,1);
+uhat = abs((2*periodogram(1:(end-1)/2+1)));
 Fs = 1/(time(2)-time(1));
 freq = Fs/2*linspace(0,1,length(u)/2+1);
 semilogy(freq(fStart:fStop)/f,real(uhat(fStart:fStop)), 'o');

SignalGeneration.m

@@ -3,10 +3,10 @@
 % "Entwicklung eines Spektrometers zur Analyse superparamagnetischer Eisenoxid-Nanopartikel für Magnetic-Particle-Imaging"
 % Gael Bringout - Oct 2014
 
-%close all
+close all
 clear all
 
-addpath('coreFunctions\')
+addpath(genpath(fullfile('.')))
 
 
 mu0         = 4*pi*1e-7; % permeability of the air.
@@ -17,17 +17,17 @@
 Ms          = 0.6/mu0; %[A/m] Magnetisation at saturation from Magnetit. From the thesis. P54
 temperature = 310;
 %signal generation and procesing
-Fs          = 25*10^6;%2.5e7; % [Hz] sampling frequency
+Fs          = 10*10^6;%2.5e7; % [Hz] sampling frequency
 fStart      = 2; % [-] staring frequency components to be displayed
-fStop       = 500; % [-] Last frequency components to be displayed. Please respect the shanon theorem :)
+fStop       = 800; % [-] Last frequency components to be displayed. Please respect the shanon theorem :)
 dt          = 1/(50*10^6); %delta used to calculate the derivative. Used to compensate the "small" sampling frequencies
 
 % MPI signal generation
 f           = 25000; % [Hz] Frequency of the signal
-nbrPeriod   = 10; % nbr of signal period
+nbrPeriod   = 4; % nbr of signal period
 time        = 0:1/Fs:nbrPeriod/f-1/Fs; % [s] time vector. Remove the last point as we start at zero.
 N           = length(time); %Number of time points
-H0          = 31831;% [A/m] Drive field amplitude (31831 A/m to generate 40 mT)
+H0          = 31831/2;% [A/m] Drive field amplitude (31831 A/m to generate 40 mT)
 B0          = H0*mu0; % [T] - H = B/mu
 
 % MPI scanner properties
@@ -39,14 +39,14 @@
 % Noise model
 kB          = 1.380650424e-23; % [J/K] Boltzmann constant 
 deltaF      = 10*10^6; % [Hz] according to Weizenecker 2007 - A simulation study...
-Rp          = 185 *10^-3; % [Ohm]  according to Weizenecker 2007 - A simulation study...
+Rp          = 0.0001*185 *10^-3; % [Ohm]  according to Weizenecker 2007 - A simulation study...
 Scaling     = 1000; % To easily scale the maximal voltage induced by the noise
 %% Particle Magnetization in function of the applied magnetic field
 
 nbrPoint = 1000;
 a=zeros(nbrPoint,1);
 mu0 = 4*pi*1e-7; % permeability of the air.
-Hpart=linspace(-10000,10000,nbrPoint);
+Hpart=linspace(-31831,31831,nbrPoint);
 Bpart = Hpart.*mu0;
 Mpart=zeros(nbrPoint,1);
 for i=1:nbrPoint,
@@ -57,7 +57,8 @@
 subplot(2,3,1)
 hold off
 plot(Bpart, Mpart, 'k');
-xlim([-max(Bpart) max(Bpart)])
+%plot(a, Mpart, 'k');
+%xlim([-max(Bpart) max(Bpart)])
 xlabel('B / (T)')
 ylabel('M / (A/m)')
 title('Particle magnetisation curve');
@@ -104,10 +105,27 @@
 %% Spectrum
 
 subplot(2,3,6)
-periodogram = fft(u);
-uhat = abs(2*periodogram(1:(end-1)/2+1));
+periodogram = fft(u)/size(u,1);
+uhat = abs((2*periodogram(1:(end-1)/2+1)));
 Fs = 1/(time(2)-time(1));
 freq = Fs/2*linspace(0,1,length(u)/2+1);
-semilogy(freq(fStart:fStop)/f,real(uhat(fStart:fStop)), 'o');
+plot(freq(fStart:fStop)/f,real(uhat(fStart:fStop)), 'o');
 title('Spectrum - # of harmonic');
-xlim([ 0 50])
+xlim([ 0 50])
+
+% filter the results
+%figure; semilogy(abs(periodogram));
+i1 = find(freq == f);
+i2 = find(freq == 2*f);
+i3 = find(freq == 3*f);
+
+AS = abs(fft(u));
+AS(i1) = AS(i1)/400;
+AS(size(u,1)-i1+2) = AS(i1);
+AS(i2) = AS(i2)/200;
+AS(size(u,1)-i2+2) = AS(i2);
+AS(i3) = AS(i3);
+AS(size(u,1)-i3+2) = AS(i3);
+
+%ybis = ifft(AS);
+%figure;plot(time,real(ybis));