Adjusting Parameters in Transcranial Doppler Audio Processing: Seeking Guidance on Gain, SWEEP, and Audio Playback

Hello there.

I'm working with an audio file, which is the output from Transcranial Doppler. The idea is to accomplish the following:

Here is the file data:

'###block_size::64

'###block_numbers::24

'###offset_binary::1984

'###session_date::20-10-2019

'###session_time::21-32-22

'###time::12-31-49-26

'###gates::2

'###prf::3003

'###samples::15726344

'###transmitting_freq1::2 MHz

'###label1::MCA_L

'###spec_gate_sv1::8

'###spec_gate_depth1::54

'###spec_gate_gain1::44

'###spec_gate_filter1::100

'###transmitting_freq2::2 MHz

Here is my code:

clc
clearvars
%% GETTING THE FLOW VELOCITY
%PRF is 3003 
% Open the binary file for reading
fid = fopen('12-31-42-34.bin');

% Read block size and block number from file header
block_size = strread(fgetl(fid),'###block_size::%d');
block_number = strread(fgetl(fid),'###block_numbers::%d');

% Skip 8 lines of header information
for k=1:8
  fgetl(fid)                                        
end

% Read the pulse repetition frequency (PRF) from the file header
prf = strread(fgetl(fid),'###prf::%d');

% Calculate the byte offset for the desired data block
offset = block_size * block_number;
fseek(fid,offset,'bof');

% Read the data from the binary file into a matrix 'a'
a = fread(fid, [4,inf],'float32')';

% Set the desire duration
duration = 5;

% Number of points for the desire duration
num_points = duration * prf;

% Set the duration of a (5 sec)
% Sample rate for the data (PRF)
Fs=prf

% Number of points for the FFT (Fast Fourier Transform)
NFFT = 128;

% Individual channels
% Get the flow velocity for each channel (assuming data format is REAL, IM, RE, IM...)
b1 = complex(a(1:num_points,2), a(1:num_points,1));  % Channel 1
b2 = complex(a(1:num_points,4), a(1:num_points,3));  % Channel 2


% Frequency vector representing the time axis
f_vec = [-floor(NFFT/2) : ceil(NFFT/2)-1] * Fs/NFFT;

% Plots
spectrogram(b1,hamming(128),64,f_vec,Fs,'yaxis');
caxis([0 20])

And I'm getting this:

I would like to replicate the same color scheme and parameters as seen in the previous image. Upon closer inspection, it appears that the gain is set to 0dB instead of 44dB, and the SWEEP is 250 fft/s. How can I incorporate these settings into my code? Additionally, I need to be able to listen to the audio output.

Thank you.

Patrícia

Answer 1

[William Rose]

@Patrícia,

For colors, you can specify a colormap that gives colors more like what you want. See here, and scroll down to see all the buil;t-in colormaps. Right now you are using parula, the default. Maybe colormap winter; will be more like what you want.

As for gain seting: Don't change the gain. The original screenshot you provided does not show units. The peaks in the original are around y=+40, and your peaks are at about +1 kHz. I suspect the y-axis units on the original screenshot are are cm/s, not kHz. A systolic velocity of 40 cm/s in the MCA is a bit below average, but not super unusual, in my experience doing TCD on adults. Just convert from kHz to cm/s.

[William Rose]

@Patrícia,

[Edit: changed the conversion equation because I left out a factor of 2.]

the conversion is

where v is positive if the target is moving toward the source, and theta is the angle between the velocity vector and the direction of illumination. We'll assume theta=0 below. Therefore the frequency shift, , is

and therefore

Let's assume c=1540 m/s (standard assumption in ultrasound scanners). Your metadata indicates f_source=2 MHz. After doing units conversions, I get

Then 1 kHz would correspond to 38.5 cm/s. Which is very similar to the peak velocity on the first image you shared above. Check my assumptions.

[Patrícia]

Thank you so much.

These are the units; that's why I thought the main problem was the gain. I'll convert to cm/s and I'll let you know.

[Patrícia]

Which variable should I apply the formula to?

[William Rose]

I would try applying the formula to b1: v1=77*b1; spectrogram(v1,hamming(128),…); But I don’t fully understand what the data in a is. You assume the columns of an are real, imaginary, which suggests the Fourier transform has already been done. If that is true, then why would you do another Fourier transform with spectrogram? Also, you said the columns of a are real, imag, but you use the opposite order when you do b1=complex(a(:,2),a(:,1)); And likewise for b2. So there’s stuff I don’t yet understand. I am using a cell phone so cannot look at this in detail. And I did not answer your question about audio output. Let me think…

[William Rose]

@Patrícia,

I suspect that my previous recommendation of how to apply the conversion formula was incorrect. I'm not certain, since I still do not understand why your data is complex. Before, I suggested that you multpliy b by the conversion factor. I now think that you should simply change the frequency axis that you use to plot the spectrogram. I would do this by getting outputs from spectrogram:

[s,f,t]=spectrogram(...);

Then I would use surf() to plot the spectrogram. And when I use surf(), instead of using the vector f (which I got from spectrogram()) for the y-axis, I would use v, where

v=77*f/1000;

for the y-axis. The factor of 1000 is there because I suspect the f vector returned by spectrogram() will have units of Hz, not kHz.

You will have to do careful testing to see if my new recommendation is correct.

For audio playback, I would try the following:

Fs=3003;        % sampling rate (Hz) 
bmag=abs(b);    % magnitude of (complex) b
sound(bmag,Fs); % play the waveform on the default sound device

Some systems won't allow a sampling rate of 3003 Hz. If so, I would resample abs(b) to a rate that is accepted by my system. So, if 22.05 kHz works on my system, I would do

newRate=22050; % audio playback sampling frequency (Hz)
bmag=resample(abs(b),newRate,3003);
sound(bmag,newRate);

Since I don't undertsand why b is complex, my sound playback recommendation may be incorrect.

I cannot check my recommendations now, since I can't run Matlab, since the license server is not responding.

[Patrícia]

The exported file from the software is a binary file (variable 'a'). For each channel, there are 2 columns, each containing a real and imaginary part. In this case, we have 2 channels, so the variable 'a' has 4 columns.

Example for Channel 1:

But I will try with this.

Thank you.

[William Rose]

@patricia, you’re welcome.

[William Rose]

@Patrícia, I'm embarrassed to say I made a mistake, which I have now corrected, in my Doppler conversion formula above. The right formula is

Sorry about that.

[William Rose]

@Patrícia, if you would post a sample file with the complex data, I would appreciate it, for learning. My TCD machine generates separate files for the positive and negative parts of velocity, but the complex vector you have is more elegant. The real and imaginary parts of your signal represent the in-phase and quadrature components of the reflected wave. As we both know, the Doppler velocity signal we see on screen includes simultaneous unequal positive and negative components. The short-time Fourier transform of a complex vector, but not of a real vector, can produce such a velocity spectrum.

I know your raw data file is big, maybe bigger than the 5 MB post-able limit, and your raw data file may contain private information. So if you just save complex vector b1 as a .mat file, and post that... Thank you.

[Patrícia]

@William Rose Don't worry. I'm using this:

c = 1540; % speed in soft tissue m/s
f_source = 2e6; %2MHz
velocity_cm_s = ((c * f*1000) / (2 * f_source))*100; 

For 1kHz, it gives me 38cm/s.

Yes, sure. So the header information that you might need is:

###block_size::64                                             
###block_numbers::24                                          
###offset_binary::1984                                                          
###gates::2                                           
###prf::3003                                            
###samples::15726344                                          
###transmitting_freq1::2 MHz                                  
###label1::MCA_L                                              
###spec_gate_sv1::8                                           
###spec_gate_depth1::54                                       
###spec_gate_gain1::44                                        
###spec_gate_filter1::100                                     
###transmitting_freq2::2 MHz                                  
###label2::MCA_R                                              
###spec_gate_sv2::8                                           
###spec_gate_depth2::57                                       
###spec_gate_gain2::38                                        
###spec_gate_filter2::100    

I added the file for the first 5 seconds. Right now, I'm trying to obtain the graph of velocity per second instead of having it in Hz.

[William Rose]

@Patrícia, Thank you for posting the file with 5 seconds of data. You wrote "I'm trying to obtain the graph of velocity per second instead of having it in Hz."

See attached script. It uses the approach I suggested above: Have spectrogram() return the power spectrum as a matrix (ps), and a vector of frequencies (f), and a vector of times (t), without making a plot. Convert the frequency vector (f) to a velocity vector (v). Then plot the power spectrum versus time and velocity, using surf().

The script makes two plots: one with frequency on the vertical axis, and one with velocity on the vertical axis (shown below), so that you can compare the results.

[Patrícia]

I was still trying to figure out how I could get there, but thank you very much! By the way, do you know any way to have the chart look more like this?

Thank you so much!

[William Rose]

@Patrícia,

Adding the command

colormap winter 

after surf(...)

gives the plot below. Not quite what you want. Look at the images of all the built-in colormaps, here, to see if anything else looks better to you.

[William Rose]

@Patrícia, if the built-in colormaps don't suit you, you can create your own, following the instructions here.

[William Rose]

@Patrícia,

For example, this shows how you can make a custom colormap with 256 colors.

patricia1=zeros(256,3);
patricia1(65:256,2)=linspace(0,1,192);
patricia1(1:65,3)=linspace(0,1,65);
patricia1(65:256,3)=linspace(1,.5,192);
save('patricia1','patricia1');

Try the new colormap.

load('patricia1');
load('patData');
figure
surf(t,v,10*log10(ps),'EdgeColor','none'); 
colormap(patricia1)
view(0,90); cb=colorbar; clim([0 20]); ylabel(cb,'Power (dB/Hz)')
xlabel('Time (s)'); ylabel('Velocity (cm/s)'); 
title('Blood Velocity')

and you get a plot using custom colormap patricia1.

If that is not what you want, then try again.

patricia2=zeros(256,3);
patricia2(128:256,2)=linspace(0,1,129);
patricia2(1:128,3)=linspace(0,1,128);
patricia2(128:256,3)=linspace(1,.5,129);
save('patricia2','patricia2');

Try it.

load('patricia2');
figure
surf(t,v,10*log10(ps),'EdgeColor','none'); 
colormap(patricia2)
view(0,90); cb=colorbar; clim([0 20]); ylabel(cb,'Power (dB/Hz)')
xlabel('Time (s)'); ylabel('Velocity (cm/s)'); 
title('Blood Velocity')

The two colormaps above go linearly from black [0,0,0] to blue [0,0,1], and then linearly from blue to SpringGreen [0,1,0.5]. Black is at row 1, SpringGreen at row 256 in both maps. By adjusting where Blue is positioned, you get more or fewer near-black colors in the image. Blue is at row 65 and 128 in patricia1 and patricia2, respectively.

[Patrícia]

Thank you for the help. If I adjust the signal gain, won't I be able to visualize it better? Because upon analyzing the first image, it seems that the signal gain was 0dB. If so, can you explain to me how I could do that?

[William Rose]

@Patrícia, the trans-cranial Doppler device has a gain setting, which you can adjust during the data collection. The manual will tell you how to change it. If you set the gain too high, the signal will become too noisy, and you will see colors, representing velocities, in the signal, where you should not see any. If the gain is too low, you will see too little color, anywhere. You will need to experiment with the gain setting to find a good value for each subject. As far as I know, the gain cannot be adjusted after the data have been collected.

[Patrícia]

The data has already been obtained, and in this case, the gain is 44dB for one of the channels, but in the software of the first image, it shows 0 dB, which suggests that they changed the gain to visualize better. So, I was thinking about how I could do that in MATLAB.

[William Rose]

@Patrícia,

There is no way to do a dB adjustment that will affect image quality after the data is collected, as far as I know. You can play with the color map, and that may make a difference. See my comments above about custom color maps, and see the matlab help for colormap. There is also a colormap editor. See the help here, and see the section "Adjust colormap to improve image detail".

[Patrícia]

Ok, thank you very much for the help. :)

[William Rose]

@Patrícia, you're welcome.