k-Wave

1. Vali
   

   Member
   Joined: Apr '20
   Posts: 3

   Hi,

   I've defined two linear transducers (2D) for plane-wave propagation.

   The first one is defined in the horizontal and the second one in vertical direction.

   If I transmit a plane-wave with the horizontal transducer and receive the signals with both, then I have a large frequency-shift in the spectrum.

   Input-Signal: Gaussian-envelope, 3 cycles and centered at 5.2MHz
   (max supported frequency of the grid ~8.3MHz)

   Spectrum of horizontal transducer: Gaussian-curve centered at ~5.2MHz
   Spectrum of vertical transducer: Gausian-curve shifted to 1.5 - 2.5MHz (channel dependent)

   Where does this large shift come from?

   Any suggestion?

   Thanks,
   Valentin

   Posted 3 years ago #

2. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

   HI Valentin,

   Sounds puzzling. Can you post your code?

   Brad.

   Posted 3 years ago #

3. Vali
   

   Member
   Joined: Apr '20
   Posts: 3

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %frequency-shift
   %
   %description:
   %Transmitting a plane wave from the horizontal elements and receive it
   %with the vertical elements (RX2). Where does this channel dependent shift
   %come from? (In compare with the spectrum of the horizontal channels,
   %where no shift occurs.)
   %note: With an increasing number of channels e.g. 192 the frequency-shift
   %will also increase.
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   clear all;
   close all;
   
   %%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %ground-settings
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   GPU = false;                                            %if GPU is available
   
   meanSpeedOfSound = 1540;                                %[m/s]
   transducerFreq = 5.2083;                                %[MHz]
   wavelength = meanSpeedOfSound/(transducerFreq*1e6);     %[m]
   
   %%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %define kgrid
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   PML_size = 20;
   PMLInside = 'false';                                    %PML is set to 'outside the grid'
   
   dxy = (wavelength)/4;                                   
   
   Ny = 243 - 2 * PML_size;                                %checkFactors() ... for computational-time
   Nx = 243 - 2 * PML_size;           
   
   Nx = 243 - 2 * PML_size;
   Ny = 243 - 2 * PML_size;
   
   kgrid = kWaveGrid(Nx, dxy, Ny, dxy);                    %make grid of given size
   
   %%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %define medium with a scatter
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   c_tissue = 1540;                                        %[m/s]
   density_tissue = 1000;                                  %[kg/m^3]
   medium.sound_speed_ref = 1540;
   medium.alpha_mode = 'no_dispersion';
   
   medium.sound_speed = c_tissue * ones(Nx, Ny);
   medium.alpha_coeff = 0.54 * ones(Nx, Ny);               %[dB/(MHz^y cm)] ... soft tissue
   medium.alpha_power = 1.5;                               %= y
   
   medium.density = density_tissue * ones(Nx, Ny);         %[kg/m^3]
   
   inclusion.dimension = [Nx/2 Ny/2 10];
   disc_c1 = round([inclusion.dimension(1,1) inclusion.dimension(1,2) inclusion.dimension(1,3)]);
   disc1 = logical(makeDisc(Nx, Ny, disc_c1(1), disc_c1(2), disc_c1(3)));
   numbOfPix1 = length(medium.sound_speed(disc1));
   lowerLimDensity = density_tissue / 1.25;
   upperLimDensity = density_tissue * 1.25;      
   
   randDensity = zeros(numbOfPix1,1);
   randDensity = (upperLimDensity - lowerLimDensity).*rand(numbOfPix1,1) + lowerLimDensity;
   medium.density(disc1) = randDensity;
   
   lowerLimSos = c_tissue / 1.25;
   upperLimSos = c_tissue * 1.25;             
   
   randSos = zeros(numbOfPix1,1);
   randSos = (upperLimSos - lowerLimSos).*rand(numbOfPix1,1) + lowerLimSos;
   medium.sound_speed(disc1) = randSos;
   
   figure(101);
   subplot(1,4,1);
   imagesc(medium.density);
   title('medium-density');
   subplot(1,4,2);
   imagesc(medium.sound_speed);
   title('medium-SoS');
   
   %%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %define input-signal
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   height = dxy*Nx;
   max_time = 2*height/max(medium.sound_speed(:));
   kgrid.t_array = makeTime(kgrid, medium.sound_speed, 0.3, max_time);
   sample_freq = 1/kgrid.dt;
   signal_freq = transducerFreq*1e6;
   num_cycles = 3;
   toneBurst_envelope = 'Gaussian';
   bandwidth = 80;
   
   %%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %define source (16 elements) and sensor (16 elements each):
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   halfDis = round(Ny/2)-8;
   source.p_mask = zeros(Nx, Ny);
   source.p_mask(20,halfDis:halfDis+15) = 1;
   sensor.mask = zeros(Nx, Ny);
   %1st RX
   sensor.mask(20,halfDis:halfDis+15) = 1;
   %1nd RX
   sensor.mask(30:45,halfDis+30) = 1;
   
   subplot(1,4,3);
   imagesc(source.p_mask);
   title('source-mask');
   subplot(1,4,4);
   imagesc(sensor.mask);
   title('sensor-mask');
   truesize;
   
   %%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %tone-burst-offset
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   phi = 0;
   toneBurst_offset = zeros(nnz(source.p_mask(:)),1);      %time-delay for elements
   toneBurst_offset = sind(phi)*[0:numel(toneBurst_offset)-1]'*dxy/meanSpeedOfSound/kgrid.dt;
   input_signal = toneBurst(sample_freq, signal_freq, num_cycles,'SignalOffset', toneBurst_offset, 'Envelope', toneBurst_envelope);
   source.p = input_signal;
   
   %%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %run simulation
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   input_args = {'PMLInside', false, 'PlotPML', false, 'DisplayMask', source.p_mask, 'PlotScale', [-0.25, 0.25]};
   if GPU, input_args = {input_args{:},'DataCast','gpuArray-single'};
   end
   sensor_data1 = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});
   if GPU, sensor_data1 = gather(sensor_data1);
   end
   
   %%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %channel selection, frequency-shift at RX2 (vertical elements)
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   ch1_RX1 = sensor_data1(1,:);
   bandwidth = 80;
   filteredSigRX1 = gaussianFilter(ch1_RX1, sample_freq, signal_freq, bandwidth, true);
   title('channel 1, RX1-horizontal');
   
   ch1_RX2 = sensor_data1(32,:);
   bandwidth = 80;
   filteredSigRX2 = gaussianFilter(ch1_RX2, sample_freq, signal_freq, bandwidth, true);
   title('channel 1, RX2-vertical');
   
   ch8_RX2 = sensor_data1(28,:);
   bandwidth = 80;
   filteredSigRX2 = gaussianFilter(ch8_RX2, sample_freq, signal_freq, bandwidth, true);
   title('channel 8, RX2-vertical');
   
   ch12_RX2 = sensor_data1(24,:);
   bandwidth = 80;
   filteredSigRX2 = gaussianFilter(ch12_RX2, sample_freq, signal_freq, bandwidth, true);
   title('channel 12, RX2-vertical');
   
   ch16_RX2 = sensor_data1(17,:);
   bandwidth = 80;
   filteredSigRX2 = gaussianFilter(ch16_RX2, sample_freq, signal_freq, bandwidth, true);
   title('channel 16, RX2-vertical');
   
   %%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %FFT-spectrum, mean value of RX1/RX2
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   testChannel = mean(sensor_data1(1:16,:));
   nfft = length(testChannel);
   ff = fft(testChannel,nfft);
   fff = ff(1:nfft/2);                      
   
   fs = 1/kgrid.dt;
   fff=fff/max(fff);
   fn = fs/2;
   df=fs/nfft;
   xfft = 0:df:fn-1;
   figure(1001);
   subplot(1,2,1);
   plot(xfft,abs(fff));
   xlabel('frequency/[Hz]');
   title('FFT mean value of RX1, normalized to 1');
   
   testChannel = mean(sensor_data1(17:32,:));
   nfft = length(testChannel);
   ff = fft(testChannel,nfft);
   fff = ff(1:nfft/2);
   fs = 1/kgrid.dt;
   fff=fff/max(fff);
   fn = fs/2;
   df=fs/nfft;
   xfft = 0:df:fn-1;                       
   
   subplot(1,2,2);
   plot(xfft,abs(fff));
   xlabel('frequency/[Hz]');
   title('FFT mean value of RX2, normalized to 1');

   Posted 3 years ago #

4. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

   You should trim the transmitted signal from the sensor data so it just contains the received signals. For example, set:

   % remove the transmitted signal
   sensor_data_trim = sensor_data1(:, 400:end);

   Then if you plot the two spectra, they cover essentially the same bandwidth. For example:

   % amplitude spectrum of average signal data
   testChannel1 = mean(sensor_data_trim(1:16,:));
   testChannel2 = mean(sensor_data_trim(17:32,:));
   [f, as1] = spect(testChannel1, 1/kgrid.dt);
   [~, as2] = spect(testChannel2, 1/kgrid.dt);
   
   % plot
   figure;
   plot(f * 1e-6, as1);
   hold on;
   plot(f * 1e-6, as2);
   legend('Sensor 1', 'Sensor 2');

   Sensor 2 has a lower amplitude because of the directional response of the finite sized sensor which you're averaging over.

   Hope that helps,

   Brad

   Posted 3 years ago #

5. Vali
   

   Member
   Joined: Apr '20
   Posts: 3

   Thanks so much for your time!

   You have totally right, I had forgotten to cut off the input-signals.
   In the previous example the frequency spectrums of both sides are equal now
   and centered around the main (center) frequency.
   But if I increase the total number of elements for Tx and Rx (e.g. 128) both spectrums
   are getting a downshift of approx. 1 MHz in center (from 5.2 to 4.2 MHz).

   Does this phenomen occur in case of frequency dependent attenuation with the depth?
   A correct TGC should compensate the frequency shift or?
   (I have tried it with the TGC settings from the B-mode example, but there is no much difference in-between.)

   >Sensor 2 has a lower amplitude because of the directional response of the finite sized sensor which you're averaging over.

   --> The sensor elements do have the same size in horizontal and vertical direction. By default the sensor elements have a
   omnidirectional sensitivity. I don't understand the difference between the amplitudes. I have also tried to align the
   sensitivity direction.

   Would you be so kind to explain that more in detail?

   Best,
   Valentin

   Posted 3 years ago #

RSS feed for this topic

Reply

You must log in to post.