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

6 次查看（过去 30 天）

显示更早的评论

Patrícia 2024-1-19

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2072031-adjusting-parameters-in-transcranial-doppler-audio-processing-seeking-guidance-on-gain-sweep-and

评论： William Rose 2024-2-16

采纳的回答： William Rose

在 MATLAB Online 中打开

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

William Rose 2024-1-19

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2072031-adjusting-parameters-in-transcranial-doppler-audio-processing-seeking-guidance-on-gain-sweep-and#answer_1393111

@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.

21 个评论
显示 19更早的评论隐藏 19更早的评论

William Rose 2024-1-20

在 MATLAB Online 中打开

@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.

William Rose 2024-1-27

在 MATLAB Online 中打开

patData.mat

@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.