%% 载入数据
clear;
clc;
close all;
% slope_rawdata = imread('-10-0-excellent-latest.bmp');lateral_location = 800; %读进来的slope为实际深度-1,lateral_location为选定的横向x位置
slope_rawdata = imread("D:\Learning Data\图像处理\第二次任务-粘弹性代码\1.bmp");
lateral_location=550;
poisition_left = 100;
poisition_right = 700;
slope_rawdata = double(slope_rawdata); %double化,范围为0-255,可用im2double,范围为0-1,两种double出的图稍微不同
slope_rawdata1 = slope_rawdata(1 : 1000,:,:);
figure(1);
imagesc(slope_rawdata); colormap(jet);
title('Raw data of slope');
colorbar;
xlabel('Lateral position (mm)');
ylabel('Time(ms)');
%% 变量定义和变量提前
multiple_t = 12; %时域倍数 频率轴间隔为1/(Aline_time*M _number),波数轴间隔为1/scan_length,为了保证两轴精度一样选择合适的倍数
Aline_time = 10 * 10^(-6); %时域采样速度
Aline_number = 750; %横向像素点
Aline_point =1024; %轴向像素点
M_number = 1000; %时域采样次数
scan_length = (0.1963 * 10 + 1.1027) * 10^(-3); %扫描范围
% multiple_x = round(2 * Aline_time * M_number * multiple_t / scan_length); %使两轴间隔相等
multiple_x = 12;
M_number_fit = M_number * multiple_t; %拟合增加采样时间
Aline_number_fit = Aline_number * multiple_x; %增加x扫描线数
scan_length_fit = scan_length * multiple_x; %增加的扫描范围,要保证x的采样率不变
lamda = 1.310; %激光器的中心波长,微米
ft = 1 / Aline_time; %时间采样频率
fx = Aline_number_fit / scan_length_fit; %空间采样频率
pre_frespec_window_f = 25 * multiple_t; %预选信号频谱窗口的频率范围
pre_frespec_window_k = 25 * multiple_x; %预选信号频谱窗口的波数范围
max_interval_left = 100; %左边最大间隔
max_interval_right = 100; %右边最大间隔
amplitude_threshold = 0.20; %方向振动最大位移的阈值
threshold_interval_left = 20; %主波左间隔阈值
threshold_interval_right = 20; %主波右间隔阈值
threshold_frespec_right = 3; %信号频谱右阈值,防止非信号频谱的干扰
threshold_frespec_down = 3; %信号频谱下阈值,防止非信号频谱的干扰
%% 数据图像优化
slope_rawdata_fft2 = fftshift(fft2(slope_rawdata1));
flog = log(abs(slope_rawdata_fft2));
filter = (flog >.53 * max(flog(:)))|(flog <.10 * max(flog(:))); %第一个数字决定修复强度,第二个数字决定弱化强度;频域里第一个数字增强信号,第二个数字抑制背景
slope_fixeddata = abs(ifft2(slope_rawdata_fft2 .* filter)); %修复之后观看修复后的slope_fixeddata有没有失真或者修复程度够不够,必看
% slope_fixeddata = slope_rawdata;
figure(2); %修复后的slope数据
imagesc(slope_fixeddata);colormap(jet);
title('fixed data of slope');
colorbar; % 添加颜色条
xlabel('Lateral position (mm)');
ylabel('Time(ms)');
%% 将slope转化为振动位移信息
mean_slope_fixed = mean(mean(slope_fixeddata(1 : 99,:)));
slope_normal = slope_fixeddata - mean_slope_fixed;
slope_amplitude_500 = slope_normal * lamda / (4 * 128);
figure(3);
imagesc(slope_amplitude_500);colormap(jet);clim([-0.20,0.20]);
colorbar; % 添加颜色条
title('Vibration amplitude');
xlabel('Lateral position (mm)');
ylabel('Time(ms)');
%% 提取局部信息
region(M_number,Aline_number) = 0;
region(:,poisition_left :poisition_right ) = 1;
slope_amplitude_500 = slope_amplitude_500 .* region;
figure(4);
imagesc(slope_amplitude_500);colormap(jet);
%% 时域增加采样点
% 时域和频谱范围
t = 0 : 1 / ft : (1 / ft * (M_number - 1)); %时间轴
t_fit = 0 : 1 / ft : (1 / ft * (M_number_fit - 1)); %增加采样时间,以提高频域精度
x = 0 : 1 / fx : (1 / fx * (Aline_number_fit - 1)); %距离轴
f = (0 : length(t_fit) / 2 - 1) * ft / length(t_fit); %频率轴
k = (0 : length(x) / 2 - 1) * fx / length(x); %波数轴
f_max = max(f); %最大频率
k_max = max(k); %最大波数
% 填充时域信息
slope_amplitude_fit(M_number_fit,Aline_number_fit) = 0;
slope_amplitude_fit(1 : M_number,1 : Aline_number) = slope_amplitude_500;
figure(4);
imagesc(slope_amplitude_fit);colormap(jet);
%% 滤波
% 二维傅里叶变换
slope_amplitude_fft2 = fftshift(fft2(slope_amplitude_fit));
Gauss = fspecial('gaussian',20,10);
[M,N] = size(Gauss);
full_Gauss(M_number_fit,Aline_number_fit) = 0;
for u = 1 : M
for v = 1 : N
full_Gauss(M_number_fit / 2 - fix(M / 2) + u,Aline_number / 2 - fix(N / 2) + v) = Gauss(u,v);
end
end
Gauss_fft2 = fftshift(fft2(full_Gauss));
fLog_filter = slope_amplitude_fft2 .* Gauss_fft2;
%fLog_filter = slope_amplitude_fft2;
fLog_filter_Log = log(abs(fLog_filter) + 1);
% 高斯滤波
%FLog_filter = imgaussfilt(fLog_filter_Log,20); %使用高斯滤波,如若不适用注释这行,取消下行注释
FLog_filter = fLog_filter_Log;
figure(5);
imagesc(FLog_filter);colormap(jet);
xlabel('Wave number (k·m^-1)'); % 横轴标签
ylabel('Frequency (kHz)'); % 纵轴标签
%% 选定信号频谱区域
FLog_leftdown = FLog_filter(M_number_fit / 2 + 1:M_number_fit / 2 + pre_frespec_window_f,Aline_number_fit / 2 + 1:-1:Aline_number_fit / 2 + 1 - pre_frespec_window_k + 1);
meanFLog_leftdown = mean(mean(FLog_leftdown));
FLog_rightdown = FLog_filter(M_number_fit / 2 + 1:M_number_fit / 2 + pre_frespec_window_f,Aline_number_fit / 2 + 1:Aline_number_fit / 2 + 1 + pre_frespec_window_k - 1);
meanFLog_rightdown = mean(mean(FLog_rightdown));
if meanFLog_leftdown > meanFLog_rightdown
entitle_frespec = FLog_leftdown;
else
entitle_frespec = FLog_rightdown;
end
figure(6);
imagesc(entitle_frespec);colormap(jet);
% 求解信号频谱
for frespec_right = 1 : pre_frespec_window_k
selected_frespec_right = entitle_frespec(:,frespec_right) > threshold_frespec_right;
right = all(~selected_frespec_right);
if right == 1
num_right = frespec_right - 1;
break;
end
if (right == 0) && (frespec_right == pre_frespec_window_k)
num_right = frespec_right;
end
end
for frespec_down = 1 : pre_frespec_window_f
selected_frespec_down = entitle_frespec(frespec_down,:) > threshold_frespec_down;
down = all(~selected_frespec_down);
if down == 1
num_down = frespec_down - 1;
break;
end
if (down == 0) && (frespec_down == pre_frespec_window_f)
num_down = frespec_down;
end
end
signal_frespec = entitle_frespec(1 : num_down,1 : num_right);
figure(7);
imagesc(signal_frespec);colormap(jet);
xlabel('Wave number (k·m ^ -1)'); % 横轴标签
ylabel('Frequency (Hz)'); % 纵轴标签
%% 频散特性曲线
pre_frequency = 2 * pre_frespec_window_f / M_number_fit * f_max;
pre_wavenumber = 2 * pre_frespec_window_k / Aline_number_fit * k_max;
cutoff_frequency = pre_frequency * num_down / pre_frespec_window_f;
cutoff_wavenumber = pre_wavenumber * num_right / pre_frespec_window_k;
f_cutoff = linspace(0,cutoff_frequency,num_down); %信号频率范围
k_cutoff = linspace(0,cutoff_wavenumber,num_right); %信号波数范围
Rol(num_down) = 0;
Col(num_down) = 0;
for j = 1 : num_down
[Rol(j),Col(j)] = max(signal_frespec(j,:)); %每个w下找最大灰度点对应的k的位置
end
% W = 2 * pi * f_cutoff;
W = f_cutoff;
K = k_cutoff(Col);
Cs = W ./ K; %相速度分布
Csstart = find((Cs > 0) & (Cs < 100),1);
f_fixed = [0,f_cutoff(Csstart : end)];
Cs_fixed = [0,Cs(Csstart : end)];
f_fit = linspace(0,cutoff_frequency,10000); %拟合频率
%p1 = polyfit(f_fixed,Cs_fixed,2); %二阶拟合
%Cs_fit1 = polyval(p1,f_fit);
Cs_fit2 = pchip(f_fixed,Cs_fixed,f_fit); %保形拟合
%Cs_fit3 = spline(f_fixed,Cs_fixed,f_fit); %样条插值
figure(8);
scatter(f_cutoff(Csstart:end),Cs(Csstart:end),'red');hold on;
%plot(f_fit,Cs_fit1);hold on;
plot(f_fit,Cs_fit2,'blue');hold on;
%plot(f_fit,Cs_fit3);hold on;
%legend('scatter','polyfit2','pchip','spline');
% 找到频散曲线比较平缓的一段
% 假设平缓段的频率范围是 [f_start, f_end]
f_start = 700; % 平缓段的起始频率
f_end = 2500; % 平缓段的结束频率
% 找到对应的频率索引
f_index_start = find(f_fit >= f_start, 1, 'first');
f_index_end = find(f_fit <= f_end, 1, 'last');
% 提取平缓段的数据
Cs_flat = Cs_fit2(f_index_start:f_index_end);
% 计算平均值和标准差
Cs_mean = mean(Cs_flat);
Cs_std = std(Cs_flat);
% 输出结果
fprintf('Average Phase Velocity: %.4f m/s\n', Cs_mean);
fprintf('Standard Deviation: %.4f m/s\n', Cs_std);
%% 选定点的波频确定
% 选定点的振动曲线
selected_point = slope_amplitude_500(:,lateral_location);
figure(9); %所求点的振动曲线
plot(selected_point);title('vibration - time curve of selected point');
xlabel('time/ms');ylabel('amplitude/um');
% 振动主序列
[minrol_amplitude,mincol_amplitude] = min(selected_point);
[maxrol_amplitude,maxcol_amplitude] = max(selected_point);
if (abs(minrol_amplitude) > maxrol_amplitude) && (maxrol_amplitude < amplitude_threshold)
for m = 1 : max_interval_left
differ_left = selected_point(mincol_amplitude - m + 1) - selected_point(mincol_amplitude - m);
differ_left_next = selected_point(mincol_amplitude - m) - selected_point(mincol_amplitude - m - 1);
if (differ_left > 0) && (differ_left_next > 0) && (m > threshold_interval_left)
interval_left = m;
break;
end
end
for n = 1 : max_interval_right
differ_right = selected_point(mincol_amplitude + n - 1) - selected_point(mincol_amplitude + n);
differ_right_next = selected_point(mincol_amplitude + n) - selected_point(mincol_amplitude + n + 1);
if (differ_right > 0) && (differ_right_next > 0) && (n > threshold_interval_right)
interval_right = n;
break;
end
end
amplitude = selected_point(mincol_amplitude - m : mincol_amplitude + n);
end
if (abs(minrol_amplitude) < maxrol_amplitude) && (abs(minrol_amplitude) < amplitude_threshold) %case 2 (unessential)
for m = 1 : max_interval_left
differ_left = selected_point(maxcol_amplitude - m + 1) - selected_point(maxcol_amplitude - m);
differ_left_next = selected_point(maxcol_amplitude - m) - selected_point(maxcol_amplitude - m - 1);
if (differ_left < 0) && (differ_left_next < 0) && (m > threshold_interval_left)
interval_left = m;
break;
end
end
for n = 1 : max_interval_right
differ_right = selected_point(maxcol_amplitude + n -1) - selected_point(maxcol_amplitude + n);
differ_right_next = selected_point(maxcol_amplitude + n) - selected_point(maxcol_amplitude + n + 1);
if (differ_right < 0) && (differ_right_next < 0) && (n > threshold_interval_right)
interval_right = n;
break;
end
end
amplitude = selected_point(maxcol_amplitude - m : maxcol_amplitude + n);
end
if (mincol_amplitude > maxcol_amplitude) && (maxrol_amplitude >= amplitude_threshold) && (abs(minrol_amplitude) >= amplitude_threshold) %case 3
for m = 1 : max_interval_left
differ_left = selected_point(maxcol_amplitude - m + 1) - selected_point(maxcol_amplitude - m);
differ_left_next = selected_point(maxcol_amplitude - m) - selected_point(maxcol_amplitude - m - 1);
if (differ_left < 0) && (differ_left_next < 0) && (m > threshold_interval_left)
interval_left = m;
break;
end
end
for n = 1 : max_interval_right
differ_right = selected_point(mincol_amplitude + n - 1) - selected_point(mincol_amplitude + n);
differ_right_next = selected_point(mincol_amplitude + n) - selected_point(mincol_amplitude + n + 1);
if (differ_right > 0) && (n > threshold_interval_right)
interval_right = n;
break;
end
end
amplitude = selected_point(maxcol_amplitude - m : mincol_amplitude + n);
end
if (mincol_amplitude < maxcol_amplitude) && (maxrol_amplitude >= amplitude_threshold) && (abs(minrol_amplitude) >= amplitude_threshold) %case 4
for m = 1 : max_interval_left
differ_left = selected_point(mincol_amplitude - m + 1) - selected_point(mincol_amplitude - m);
differ_left_next = selected_point(mincol_amplitude - m) - selected_point(mincol_amplitude - m - 1);
if (differ_left > 0) && (differ_left_next > 0) && (m > threshold_interval_left)
interval_left = m;
break;
end
end
for n = 1 : max_interval_right
differ_right = selected_point(maxcol_amplitude + n - 1) - selected_point(maxcol_amplitude + n);
differ_right_next = selected_point(maxcol_amplitude + n) - selected_point(maxcol_amplitude + n + 1);
if (differ_right < 0) && (differ_right_next < 0) && (n > threshold_interval_right)
interval_right = n;
break;
end
end
amplitude = selected_point(mincol_amplitude - m : maxcol_amplitude + n);
end
figure(10);
plot(amplitude);
提取amplitude的主频
最新发布