（附超声射频数据）matlab实现b超原始信号数据生成图像

本文来自09年川大一位同学写的实验报告，算法实现很粗糙，适合新手了解B超信号处理过程。
实验数据见百度云
链接：https://pan.baidu.com/s/1Ip0Ecnm26Ol2-b7HAg_65A
提取码：0e33

1. 超声RF信号的获取与导入

在实验室中已经获得了超声的rf信号，这些数据存放在采集程序生成的.dat文件中，这些信号是超声生成图像的原始数据。
在数据处理的过程中，需要将.dat文件中的数据导入到matlab中，这需要了解.dat文件中数据的安排规则:
.dat文件的数据格式为：16bit（每个点占用位宽，有符号数） ×8035（点数）×240线（线数），共3085400字节。先从小到大（0～8034）存第1线8035点，然后存第2线8035点，以此类推，直到240线结束。整个*.Dat文件是裸数据，没有别的信息。
Matlab中有打开文件的函数为fopen，对文件中的数据进行读取的函数为fread，则将.dat文件中的数据读到变量中。需要注意的是，读入的数据的格式为一维向量，根据数据格式，将这一维向量变换为240×7618的矩阵。则可以准确得到一帧数据。

function Data = changeB2T(str)
%将二进制的数据转为十进制，并且存在数组Data中，str为每一帧图像数据的二进制文件的文件名
fid = fopen(str,'r');
Data = fread(fid,'short');
%data是一个n*1的矩阵，现按照超声生成的数据格式转换为m*n矩阵，其中m代表扫描线数目1~240，n代表每一条扫描线的数据，1~7618
Data = Data';   %转置
for i = 1:240
    da(i,:) = Data(7618*(i-1)+1:7618*i);
end
Data = da;
fclose(fid);

一线数据的时域显示

2.超声信号的预处理

由于得到的数据可能混有噪声等，所以需要对每一线的数据进行滤波、抽取处理。

2.1、数据的频移
超声发射的中心频率可以将数据的频域中找到，而我们需要保留的就是中心频率周围的信息，其他频率的数据需要滤除。
对数据进行滤波的思路为：将每一线的数据在其频域中找到中心频率，然后将中心频率移到零频附近，对移频后的数据进行低通滤波，可以得到滤波之后的数据。
移频的过程是：在时域中，将信号的实部乘以 ，虚部乘以 ，其中， 是频域中中心频率对应的点数， 是一线数据一共的数据点数。
需要计算的是从数据的频域中找到N，方法如下：由于数据在频域中心对称，所以，只需要找到频域数据从[1,n/2]个数据中最大值对应的点数，即是N。

mag = abs(fft(Data));
m = size(mag,2);
[max k] = max(mag(1:fix(m/2)));
w = -2*pi*k/N;
%%%%%%%%%%%%oiginal
x = 0:N-1;
Q = cos(w*x).*Data;
I = sin(w*x).*Data;

2.2、数据的滤波
数据进行频移之后，则有意义的数据被移至低频段，我们可以直接进行低通滤波即可。
在matlab中，一维低通滤波可以通过以下两个函数进行：fir = fir1(n,wn) 可以构造线性一维滤波器，其中n为滤波器的阶数，wn为截止频率；fftfilt(fir,Q)可以对Q中的数据用fir滤波器进行滤波。

%下面开始设计低通滤波器，采样频率为1/T ,归一化截止频率为wn，截止频率为f
%效果不太明显，还需改进,抽取之后出现明显混叠，抽取的后的采样频率为2.5MHZ，所以截止频率必须在1MHZ以内
T = 1/40000;
f = 100;
wn = 2*T*f*1.001;  
n = 70;
fir = fir1(n,wn);


%下面开始滤波
Qf = fftfilt(fir,Q);   %滤波后的时域信号,实部
If = fftfilt(fir,I);   %滤波后的时域信号，虚部
%对实部和虚部进行滤波后，将实部和虚部的数据直接进行合成，得到滤波后的数据：
Fdata = sqrt(Q.^2 + I.^2);

2.3、rf信号的抽取
由于数据量太大，需要对数据进行抽取，在时域内对数据进行抽取，则频域内数据频谱的间隔距离增大，可能会造成频谱混叠，所以，在抽取时必须保证抽取后的采样频率仍然大于我们所需要的中心频率的两倍。
我们采集数据的采样频率为40MHz，在数据的滤波过程中，设计的滤波器截止频率为100Hz，所以，我们设定抽取因子D=16，则采样频率变为40/6 = 6 MHz，可以看出，6MHz远远大于截止频率100Hz，所以，D=16的抽取因子不会造成频谱混叠。

%对数据进行抽取，每16个点抽取一个点,得到的数据在sdata里面
index = 1;
p = 16*(index-1)+1;
while p<N
   Sdata(index) = Fdata(1,p);
   p = 16*(index-1)+1;
   index = index + 1;
end

预处理总程序

function Sdata = processF(Data)
%对每一线的数据进行处理，包括滤波、对数线性变换
%Data是超声获得的每一条扫描线上的数据，数据格式为1*7618
%中心频率3.5MHZ对应的点数为K,总点数为7618，超声采样频率为40MHZ,频移的频率为w
%Sdata是时域信号
N =7618;
M = 3.5;
C = 40;
%k = fix(N*M/C);
mag = abs(fft(Data));
m = size(mag,2);
[p k] = max(mag(1:fix(m/2)));
w = -2*pi*k/N;
%%%%%%%%%%%%oiginal
x = 0:N-1;
Q = cos(w*x).*Data;
I = sin(w*x).*Data;

%下面开始设计低通滤波器，截止频率稍大于1w，采样频率为1/T ,归一化截止频率为wn，截止频率为f
%效果不太明显，还需改进,抽取之后出现明显混叠，抽取的后的采样频率为2.5MHZ，所以截止频率必须在1MHZ以内
T = 1/40000;
f = 100;
wn = 2*T*f*1.001;  
n = 70;
fir = fir1(n,wn);


%下面开始滤波
Qf = fftfilt(fir,Q);   %滤波后的时域信号,实部
If = fftfilt(fir,I);   %滤波后的时域信号，虚部
Fdata = sqrt(Q.^2 + I.^2);
%需要把频域数据中心化表示，用函数fftshift(频域数据)

%对数据进行抽取，每16个点抽取一个点,得到的数据在sdata里面
index = 1;
p = 16*(index-1)+1;
while p<N
   Sdata(index) = Fdata(1,p);
   p = 16*(index-1)+1;
   index = index + 1;
end
Sdata = Sdata/255;
%下面对滤波抽取之后的数据进行检测，结果说明数据的处理是正确的
% figure;
% plot(abs(fftshift(fft(Sdata))));
%下面对每一线的数据进行灰度对数变换
k = size(Sdata,2);
a = 1.5;
b = 0.02;

for n=1:k
    Sdata(n) = a*log(Sdata(n)+1) + b;
    if Sdata(n) > 255
        Sdata(n) = 255;
    end
end

3.图像的生成

function CData = frameProcess(str)
%对每一帧的数据进行处理，输出的是每一帧的完全处理好的数据
data = changeB2T(str);
%对数据进行绝对值处理
data = abs(data);
%对每一线的数据进行处理
for n = 1:240
    Cdata(n,:) = processF(data(n,:));
end
%进行扫描线重排
for n = 186:240
    CData(n-185,:) = Cdata(n,:);
end
for n = 56:240
    CData(n,:) = Cdata(n-55,:);
end
%转置
CData = CData';

转换为扇形的图像

主程序

%将每一帧处理好的数据进行整合

d1 = frameProcess('SongYin_LV1.dat');
d2 = frameProcess('SongYin_LV2.dat');
d3 = frameProcess('SongYin_LV3.dat');
d4 = frameProcess('SongYin_LV4.dat');
d5 = frameProcess('SongYin_LV5.dat');
% 对还未进行坐标变换的图像进行检测，看是否与标准图像类似
d = d1 + d2 + d3 + d4 + d5;
d = d/5;

%d = d1;
%d中是有效数据
imshow(d);
%下面进行坐标变换，从直角坐标系转换到扇形坐标系下
%设新的图像的大小仍为原数据矩阵大小,则每线有m个数据，一共有n线
[m n] = size(d);
data = zeros(m,n);
data = data - 1;
%根据仪器构造，扫描角度为68°，设角度圆心到图像的起始半径为r（像素），整个扇形有效数据的半径为n+r（像素）
deltasita = 68/240/180*pi;
Fsita = 68/180*pi;
r = 70;
Y = r*cos(34/180*pi);
X = fix(n/2);
for sita = deltasita:deltasita:Fsita
    for i = 1:m
       R = r + i;
       alpha = 236/180*pi + sita;      
       x = R*cos(alpha);
       y = R*sin(alpha);
       x = fix(x + X);
       y = -fix(y - Y);
       if (x<1 || x>n) continue;
       end
       if (y<1 || y>m) continue;      
       end
       if data(y,x) == -1
           data(y,x) = d(i,fix(sita/deltasita));
       end
    end
end
% %开始插值
% for i = 1:m
%     for j = 1:n
%         if data(i,j) == -1
%             x = j - X;   %扇形坐标系下坐标
%             y = - Y - i;
%             tansita = y/x;
%             sita = atan(tansita);   %角度弧度
%             if sita < 0
%                 sita = sita + 2*pi;
%             else
%                 sita = sita + pi;
%             end
%             if (sita <= 236/180*pi) || (sita>=304/180*pi)
%                 data(i,j) = 0;
%                 continue;    %不在有效扇形区域内，跳过
%             end
%             sita = sita - 236/180*pi;
%             R = fix(abs(x/cos(sita)));
%             n = fix(sita/deltasita);
%             if (n>240) || (n<1) 
%                 continue;
%             end
%             if (R >478) || (R<1)
%                 continue;
%             end
%             data(i,j) = d(R,n);
%         end
%     end
% end
imshow(data);

萌新可以看看这篇翻译文章，对b超原理讲的很透彻。
简单解释超声波成像的工作原理
https://zhuanlan.zhihu.com/p/359077343