Changeset 780
 Timestamp:
 May 28, 2014, 9:34:23 PM (10 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/series/beam_forming.m
r779 r780 12 12 ParamOut.ProjObject='off';%can use projection object(option 'off'/'on', 13 13 ParamOut.Mask='off';%can use mask option (option 'off'/'on', 'off' by default) 14 ParamOut.OutputDirExt='.p_formed ';%set the output dir extension14 ParamOut.OutputDirExt='.p_formed_1';%set the output dir extension 15 15 ParamOut.OutputFileMode='NbInput';% '=NbInput': 1 output file per input file index, '=NbInput_i': 1 file per input file index i, '=NbSlice': 1 file per slice 16 16 %check the input files … … 96 96 Data.VarDimName={'Coord_x','Coord_z',{'Coord_z','Coord_x'}}; 97 97 Data.Coord_x=5*(nbvoie_reception0.5)/numel(nbvoie_reception); % totql length of e 98 Data.Coord_z=(1:A)/133 ;% to check from input parameter .... 98 99 %%%%%% 99 100 % 100 101 % while test_fin_fichier>0 101 102 % if read_data==1 102 directory='manip_lgit';%%%%%%%%%%%%%%%%%103 name='test';%%%%%%%%%%%%%%%%%103 %directory='manip_lgit';%%%%%%%%%%%%%%%%% 104 %name='test';%%%%%%%%%%%%%%%%% 104 105 % number=2; 105 106 numero_tir_debut=1;107 numero_tir_fin=numero_tir_fin_old+pas_fichier;108 109 % eval(['load ' directory '\' name '.mat'])110 matrice_finale=zeros(A,length(nbvoie_reception),numero_tir_fin);111 time=(b/rsf+[0:A1]/rsf);112 freq1=0.5;freq2=1.5;113 [BB AA]=butter(4,[freq1 freq2]/rsf*2);114 106 number=1; %subsequence index (from 1 to 5) 115 OutputDirExt=['.p_formed_' num2str(number)]; %TODO: loop on 'number' from 1 to 5 116 for ii=1:length(nbvoie_reception)%=64 117 %ii 118 %eval(['fid=fopen(''E:\ManipLGITLecoeur\' directory '\' name '_' num2str(number) '_' num2str(nbvoie_reception(ii)) '.dat'',''r'');']); 119 filename=fullfile_uvmat(RootPath,SubDir,RootFile,FileExt,NomType,number,ii) 120 fid=fopen(filename); 121 toto=zeros(Nsequence*A*numero_tir_fin+31,1); 122 toto=fread(fid,numero_tir_fin*A*Nsequence+31,'int16','ieeele') ; 123 toto=double(bitxor(uint16(toto),uint16(2048))); 124 toto(1:31)=[];toto(numero_tir_fin*A*Nsequence)=mean(toto); 125 fclose(fid); 126 127 tata=reshape(toto2048,A,numero_tir_fin); 128 matrice_finale(:,ii,:)=reshape(tata,[A,1,numero_tir_fin]); 129 clear toto tata 130 end 131 132 matrice_finale(:,:,numero_tir_debut:numero_tir_fin_old)=[]; 133 numero_tir_fin=numero_tir_fin1; 134 matrice_finale=reshape(filtfilt(BB,AA,matrice_finale(:,:)),size(matrice_finale)); 135 136 % if soustraction==1 137 % eval(['load moyenne_' name '_' num2str(number) '.mat matrice_finale_moy']) 138 % for kk=1:size(matrice_finale,3) 139 % matrice_finale(:,:,kk)=matrice_finale(:,:,kk)matrice_finale_moy; 140 % end 141 % end 142 143 eval(['save matrice_finale_' num2str(numero_tir_fin_old) '_' num2str(numero_tir_fin) '.mat']) 144 145 %%%%%%%%%%%%%%Imagerie 146 fe=rsf*1e6; 147 cc=1475; 148 hanning_window=25; 149 hanning_vect=hanning(2*hanning_window+1); 150 interval=[1:size(matrice_finale,1)]; 151 freq=0:fe/length(interval):fe*(11/length(interval)); 152 153 pas_reseau_z=0.75e3;%0.75e3 154 pas_reseau_r=0; 155 voie_mean=length(nbvoie_reception)/2;%32; 156 reseau_z=[0:length(nbvoie_reception)1]*pas_reseau_z; 157 reseau_z=reseau_zreseau_z(voie_mean); 158 reseau_r=[0:length(nbvoie_reception)1]*pas_reseau_r; 159 reseau_r=reseau_rreseau_r(voie_mean); 160 161 debut_r=(time(1)+20)*1e6*cc/2; 162 fin_r=(time(end)20)*1e6*cc/2; 163 164 image_r=debut_r:.5e3:fin_r; 165 image_z=24e3:.75e3:24e3; 166 167 image_fin=zeros(length(image_r),length(image_z),size(matrice_finale,3)); 168 %image_fin_bis=zeros(length(image_r),length(image_z),size(matrice_finale,3)); 169 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 170 Data.Coord_z=(1:size(image_fin,1))/133 ;% to check from input parameter .... 171 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 172 for kk=1:size(matrice_finale,3) 173 kk 174 tir=kk; 175 signal=squeeze(matrice_finale(interval,:,tir)); 176 tata_fft=zeros(size(signal,1),size(signal,2)); 177 tata_fft=fft(signal,[],1); 178 179 if kk==1 107 OutputDirExt=['.p_formed_' num2str(number)]; %TODO: loop on 'number' from 1 to 5 108 numero_tir_fin_old=1; 109 pas_fichier=20;% nbre of successive shoots to read (to account for computer memory limit) 110 Nmoy=20; %%%%% value 20 FOR TEST : to shift to VALUE 8000 set by the .mat file 111 112 test_fin_fichier=1;% test to stop input file reading 113 while test_fin_fichier>0 114 numero_tir_debut=1; 115 numero_tir_fin=numero_tir_debut+pas_fichier; 116 117 % eval(['load ' directory '\' name '.mat']) 118 matrice_finale=zeros(A,length(nbvoie_reception),numero_tir_fin);%A=nbre of times (coord z)=2650, numero_tir_fin=time index 119 time=(b/rsf+[0:A1]/rsf); %b=250, rsf=10, 120 freq1=0.5;freq2=1.5; 121 [BB AA]=butter(4,[freq1 freq2]/rsf*2); 122 123 for ii=1:length(nbvoie_reception)%=64 124 %eval(['fid=fopen(''E:\ManipLGITLecoeur\' directory '\' name '_' num2str(number) '_' num2str(nbvoie_reception(ii)) '.dat'',''r'');']); 125 filename=fullfile_uvmat(RootPath,SubDir,RootFile,FileExt,NomType,number,[],ii); 126 fid=fopen(filename); 127 toto=zeros(Nsequence*A*numero_tir_fin+31,1);% Nsequence=1 128 toto=fread(fid,numero_tir_fin*A*Nsequence+31,'int16','ieeele') ; 129 toto=double(bitxor(uint16(toto),uint16(2048))); 130 toto(1:31)=[];toto(numero_tir_fin*A*Nsequence)=mean(toto); 131 fclose(fid); 180 132 181 matrice_freq_mean=mean(abs(fft(signal,[],1)),2); 182 X=[freq1*1e6 freq2*1e6]; 183 [I J]=find(freq>=X(1) & freq<=X(2)); 184 int_freq=find(matrice_freq_mean(round(1:length(freq)/2))>max(matrice_freq_mean(round(1:length(freq)/2)))/2); 185 bandwidth=freq(int_freq(end)int_freq(1)); 186 %clear matrice_freq_mean 187 end 188 189 for ii=1:length(image_r) 190 [ii kk] 191 for jj=1:length(image_z) 133 tata=reshape(toto2048,A,numero_tir_fin); 134 matrice_finale(:,ii,:)=reshape(tata,[A,1,numero_tir_fin]); 135 clear toto tata 136 end 137 138 matrice_finale(:,:,numero_tir_debut:numero_tir_fin_old)=[]; 139 numero_tir_fin=numero_tir_fin1; 140 matrice_finale=reshape(filtfilt(BB,AA,matrice_finale(:,:)),size(matrice_finale)); 141 142 % if soustraction==1 143 % eval(['load moyenne_' name '_' num2str(number) '.mat matrice_finale_moy']) 144 % for kk=1:size(matrice_finale,3) 145 % matrice_finale(:,:,kk)=matrice_finale(:,:,kk)matrice_finale_moy; 146 % end 147 % end 148 %eval(['save matrice_finale_' num2str(numero_tir_fin_old) '_' num2str(numero_tir_fin) '.mat']) 149 150 %%%%%%%%%%%%%%Imagerie 151 fe=rsf*1e6;%acquisition frequency 152 cc=1475;%speed of sound 153 hanning_window=25; 154 hanning_vect=hanning(2*hanning_window+1); 155 interval=[1:size(matrice_finale,1)]; 156 freq=0:fe/length(interval):fe*(11/length(interval)); 157 158 pas_reseau_z=0.75e3;%0.75e3 159 pas_reseau_r=0; 160 voie_mean=length(nbvoie_reception)/2;%32; 161 reseau_z=[0:length(nbvoie_reception)1]*pas_reseau_z; 162 reseau_z=reseau_zreseau_z(voie_mean); 163 reseau_r=[0:length(nbvoie_reception)1]*pas_reseau_r; 164 reseau_r=reseau_rreseau_r(voie_mean); 165 166 debut_r=(time(1)+20)*1e6*cc/2; 167 fin_r=(time(end)20)*1e6*cc/2; 168 169 image_r=debut_r:.5e3:fin_r; 170 image_z=24e3:.75e3:24e3; 171 172 image_fin=zeros(length(image_r),length(image_z),size(matrice_finale,3));%size=(332,65,pas_fichier) 173 %image_fin_bis=zeros(length(image_r),length(image_z),size(matrice_finale,3)); 174 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 175 176 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 177 for kk=1:size(matrice_finale,3) 178 kk 179 tir=kk; 180 signal=squeeze(matrice_finale(interval,:,tir)); 181 tata_fft=zeros(size(signal,1),size(signal,2)); 182 tata_fft=fft(signal,[],1); 183 184 if kk==1 192 185 193 delay=zeros(length(nbvoie_reception),1); 194 delay=1/cc*sqrt((reseau_zimage_z(jj)).^2+(reseau_rimage_r(ii)).^2); 186 matrice_freq_mean=mean(abs(fft(signal,[],1)),2); 187 X=[freq1*1e6 freq2*1e6]; 188 [I J]=find(freq>=X(1) & freq<=X(2)); 189 int_freq=find(matrice_freq_mean(round(1:length(freq)/2))>max(matrice_freq_mean(round(1:length(freq)/2)))/2); 190 bandwidth=freq(int_freq(end)int_freq(1)); 191 %clear matrice_freq_mean 192 end 193 194 for ii=1:length(image_r) 195 [ii kk] 196 for jj=1:length(image_z) 197 198 delay=zeros(length(nbvoie_reception),1); 199 delay=1/cc*sqrt((reseau_zimage_z(jj)).^2+(reseau_rimage_r(ii)).^2); 200 201 [ind centre_z]=min(abs((reseau_zimage_z(jj)))); 202 interval_utile=round(((delay(centre_z)+1/cc*abs(image_r(ii)))*fe)(b+interval(1)1)+round(length(motifbase)/2)+[fe/bandwidth/2:fe/bandwidth/2]); 203 delay=delaydelay(centre_z); 204 205 hanning_vecteur=zeros(1,length(nbvoie_reception)); 206 if centre_z>hanning_window & centre_z<(length(nbvoie_reception)hanning_window) 207 hanning_vecteur(centre_z+[hanning_window:hanning_window])=hanning_vect; 208 elseif centre_z<=hanning_window 209 test=hanning_vect((centre_z+[hanning_window:hanning_window])>=1); 210 hanning_vecteur(1:length(test))=test; 211 elseif centre_z>=(length(nbvoie_reception)hanning_window) 212 test=hanning_vect((centre_z+[hanning_window:hanning_window])<=length(nbvoie_reception)); 213 hanning_vecteur(length(nbvoie_reception)+[length(test)+1:0])=test; 214 end 215 hanning_vecteur=hanning_vecteur/norm(hanning_vecteur); 216 clear test; 217 218 amplitude_weight=ones(size(signal,1),1)*hanning_vecteur; 219 signal_new_rec=zeros(size(signal,1),length(nbvoie_reception)); 220 221 tata=zeros(size(signal,1),size(signal,2)); 222 tata(J,:)=tata_fft(J,:).*exp(1i*2*pi*(freq(J)'*delay)); 223 signal_new_rec=2*real(ifft(tata,[],1)).*amplitude_weight; 224 index_interval_utile=find(interval_utile>0 & interval_utile<size(signal,1)); 225 toto=zeros(length(index_interval_utile),1); 226 toto=mean(signal_new_rec(interval_utile(index_interval_utile),:),2); 227 image_fin(ii,jj,kk)=sqrt(mean(toto.^2)); 228 clear signal_bis interval_utile index_interval_utile hanning_vecteur 229 end 230 end 231 end 232 233 clear signal_new_em signal_new_rec m delay toto toto_bis tata tata_fft 234 235 if affichage==1 236 for kk=1:size(image_fin,3) 195 237 196 [ind centre_z]=min(abs((reseau_zimage_z(jj)))); 197 interval_utile=round(((delay(centre_z)+1/cc*abs(image_r(ii)))*fe)(b+interval(1)1)+round(length(motifbase)/2)+[fe/bandwidth/2:fe/bandwidth/2]); 198 delay=delaydelay(centre_z); 199 200 hanning_vecteur=zeros(1,length(nbvoie_reception)); 201 if centre_z>hanning_window & centre_z<(length(nbvoie_reception)hanning_window) 202 hanning_vecteur(centre_z+[hanning_window:hanning_window])=hanning_vect; 203 elseif centre_z<=hanning_window 204 test=hanning_vect((centre_z+[hanning_window:hanning_window])>=1); 205 hanning_vecteur(1:length(test))=test; 206 elseif centre_z>=(length(nbvoie_reception)hanning_window) 207 test=hanning_vect((centre_z+[hanning_window:hanning_window])<=length(nbvoie_reception)); 208 hanning_vecteur(length(nbvoie_reception)+[length(test)+1:0])=test; 209 end 210 hanning_vecteur=hanning_vecteur/norm(hanning_vecteur); 211 clear test; 212 213 amplitude_weight=ones(size(signal,1),1)*hanning_vecteur; 214 signal_new_rec=zeros(size(signal,1),length(nbvoie_reception)); 215 216 tata=zeros(size(signal,1),size(signal,2)); 217 tata(J,:)=tata_fft(J,:).*exp(1i*2*pi*(freq(J)'*delay)); 218 signal_new_rec=2*real(ifft(tata,[],1)).*amplitude_weight; 219 index_interval_utile=find(interval_utile>0 & interval_utile<size(signal,1)); 220 toto=zeros(length(index_interval_utile),1); 221 toto=mean(signal_new_rec(interval_utile(index_interval_utile),:),2); 222 image_fin(ii,jj,kk)=sqrt(mean(toto.^2)); 223 clear signal_bis interval_utile index_interval_utile hanning_vecteur 224 end 225 end 226 end 227 228 clear signal_new_em signal_new_rec m delay toto toto_bis tata tata_fft 229 230 if affichage==1 231 for kk=1:size(image_fin,3) 232 233 figure(1) 234 imagesc(image_r*1e2,image_z*1e2,image_fin(:,:,kk)'/max(max(image_fin(:,:,kk)))'); 235 title(['avec beamforming  energie max = ' num2str(max(max(image_fin(:,:,kk))))]) 236 colorbar; 237 xlabel('r (cm)');ylabel('z (cm)'); 238 drawnow 239 pause(.2); 240 end 241 end 242 243 clear matrice_finale 244 245 %%%%%%% TO ADAPT 246 for iii=1:size(image_fin,3) 247 Data.Pressure=image_fin(:,:,iii); 248 FileIndex=iii+((kk1)*nombre_tir_relu)+((qq1)*Nmoy);%%%%%%TO CHECK!!!!! 249 %%%%%%%%%% 250 %eval(['save analyse_' name '_' num2str(number) '_' num2str(numero_tir_fin_old) '_' num2str(numero_tir_fin) '.mat']) 251 OutputDir=[Param.OutputSubDir OutputDirExt];% subdirectory for output files 252 OutputFile=fullfile_uvmat(OutputDir,'','signal','.nc','_00001',FileIndex); 253 error=struct2nc(OutputFile,Data);%save result file 254 if isempty(error) 255 disp(['output file ' OutputFile ' written']) 256 else 257 disp(error) 258 end 259 end 260 %numero_tir_fin_old=numero_tir_fin+1; 261 238 figure(1) 239 imagesc(image_r*1e2,image_z*1e2,image_fin(:,:,kk)'/max(max(image_fin(:,:,kk)))'); 240 title(['avec beamforming  energie max = ' num2str(max(max(image_fin(:,:,kk))))]) 241 colorbar; 242 xlabel('r (cm)');ylabel('z (cm)'); 243 drawnow 244 pause(.2); 245 end 246 end 247 248 clear matrice_finale 249 250 %%%%%%% TO ADAPT 251 for iii=1:size(image_fin,3) 252 Data.Pressure=image_fin(:,:,iii); 253 FileIndex=numero_tir_fin  pas_fichier+iii;%%%%%%TO CHECK!!!!! 254 %%%%%%%%%% 255 %eval(['save analyse_' name '_' num2str(number) '_' num2str(numero_tir_fin_old) '_' num2str(numero_tir_fin) '.mat']) 256 OutputDir=[Param.OutputSubDir OutputDirExt];% subdirectory for output files 257 OutputFile=fullfile_uvmat(OutputDir,'','signal','.nc','_00001',FileIndex); 258 error=struct2nc(OutputFile,Data);%save result file 259 if isempty(error) 260 disp(['output file ' OutputFile ' written']) 261 else 262 disp(error) 263 end 264 end 265 numero_tir_fin_old=numero_tir_fin+1; 266 if (numero_tir_fin_old+pas_fichier)>Nmoy 267 test_fin_fichier=1; 268 end 269 end 270 271
Note: See TracChangeset
for help on using the changeset viewer.