现在已有两个函数如下:
%% 修改后的 generate_tf_waterfall 函数
function [composite_signal, t_total, Fs] = generate_tf_waterfall(grayImg, f1, f2, t1, t2)
% 参数验证
if nargin < 5
error('需要5个输入参数');
end
if f1 >= f2
error('f1 必须小于 f2');
end
if t1 >= t2
error('t1 必须小于 t2');
end
if ~ismatrix(grayImg) || ~isnumeric(grayImg)
error('grayImg 必须是二维数值数组');
end
if any(grayImg(:) > 1 | grayImg(:) < 0)
warning('二值图像应在[0,1]范围内,自动二值化处理');
grayImg = grayImg > 0.5;
end
% === 关键修改1: 自适应采样率 ===
k = 2.5; % 采样率倍数 (满足 Nyquist 定律)
Fs = max(1000, k * f2); % 确保最低采样率1000Hz且Fs > 2*f2
if Fs <= 2*f2
Fs = 2.1 * f2; % 强制满足 Nyquist
warning('调整Fs=%.1fHz以满足Nyquist要求', Fs);
end
% 获取图像尺寸
[M, N] = size(grayImg);
total_duration = t2 - t1;
% 计算时间步长
delta_t = total_duration / max(N, 1); % 避免除零
% 创建总时间向量
t_total = linspace(t1, t2, ceil(total_duration * Fs));
composite_signal = zeros(1, numel(t_total));
% 预计算频率映射
freq_bins = linspace(f2, f1, M+1); % 顶部行=高频, 底部行=低频
% 处理每个时间列
for j = 1:N
% 当前时间区间
t_start_pixel = t1 + (j-1) * delta_t;
t_end_pixel = t1 + j * delta_t;
non_zero_rows = find(grayImg(:, j));
for k = 1:numel(non_zero_rows)
row_idx = non_zero_rows(k);
% 计算频率范围
if row_idx < M
f_low = freq_bins(row_idx + 1);
f_high = freq_bins(row_idx);
else
f_low = freq_bins(end-1);
f_high = freq_bins(end);
end
% 频率边界保护
f_low = max(f1, min(f2, f_low));
f_high = max(f1, min(f2, f_high));
if f_low >= f_high
f_high = f_low + 0.1*(f2-f1); % 最小带宽保护
end
% === 关键修改2: 传递Fs参数 ===
[pixel_signal, t_pixel] = generate_tf_rectangle(...
t_start_pixel, t_end_pixel, f_low, f_high, Fs); % 传入Fs
if isempty(pixel_signal)
continue;
end
% 对齐时间向量并叠加信号
idx_start = find(t_total >= t_pixel(1), 1);
if isempty(idx_start), idx_start = 1; end
idx_end = find(t_total <= t_pixel(end), 1, 'last');
if isempty(idx_end), idx_end = numel(t_total); end
valid_idx = idx_start:min(idx_end, numel(composite_signal));
sig_len = length(pixel_signal);
valid_len = length(valid_idx);
if valid_len > 0
if sig_len > valid_len
composite_signal(valid_idx) = composite_signal(valid_idx) + ...
pixel_signal(1:valid_len);
else
composite_signal(valid_idx(1:sig_len)) = ...
composite_signal(valid_idx(1:sig_len)) + pixel_signal;
end
end
end
% 进度显示
if mod(j, 10) == 0
fprintf('处理进度: %d/%d 列 (%.1f%%)\n', j, N, j/N*100);
end
end
% 归一化信号
if ~isempty(composite_signal)
composite_signal = composite_signal / max(abs(composite_signal));
end
% 可视化
if nargout == 0 || ~isempty(composite_signal)
visualize_spectrogram(composite_signal, Fs, grayImg, f1, f2, t1, t2);
end
end
%% 修改后的 generate_tf_rectangle 函数
function [signal, t] = generate_tf_rectangle(t_start, t_end, f_low, f_high, Fs)
% 参数验证
if nargin < 5
error('需要5个输入参数');
end
if t_start >= t_end
error('t_start 必须小于 t_end');
end
if f_low >= f_high
error('f_low 必须小于 f_high');
end
% === 关键修改3: 动态Nyquist校验 ===
nyquist = Fs/2;
if f_high > nyquist
f_high = min(nyquist * 0.95, f_high); % 强制低于Nyquist频率
f_low = min(f_high - 1, f_low); % 保持最小带宽
warning('频率(%.1f-%.1fHz)超出Nyquist(%.1fHz), 已调整', ...
f_low, f_high, nyquist);
end
min_samples = 36; % filtfilt要求的最小样本数
buffer_time = 0.1; % 前后缓冲时间(秒)
% 计算原始时间区间长度
orig_duration = t_end - t_start;
orig_samples = ceil(orig_duration * Fs);
% === 关键修改4: 优化时间向量范围 ===
start_time = max(0, t_start - buffer_time);
end_time = t_end + buffer_time;
t = start_time:1/Fs:end_time;
signal = zeros(size(t));
% 找到有效时间区间
idx_start = find(t >= t_start, 1);
if isempty(idx_start), idx_start = 1; end
idx_end = find(t <= t_end, 1, 'last');
if isempty(idx_end), idx_end = length(t); end
valid_idx = idx_start:idx_end;
segment_length = length(valid_idx);
% 处理极短时间区间
if segment_length < min_samples
if segment_length > 10
center_freq = (f_low + f_high)/2;
[signal_seg, t_seg] = generate_single_tone(t_start, t_end, center_freq, Fs);
else
[signal_seg, t_seg] = generate_impulse((t_start+t_end)/2, f_low, f_high, Fs);
end
% 将短信号放入对应区间
[~, pos] = min(abs(t - t_seg(1)));
valid_pos = pos:min(pos+length(signal_seg)-1, length(signal));
signal(valid_pos) = signal_seg(1:length(valid_pos));
return;
end
% 生成高斯白噪声
noise = randn(1, segment_length);
% 设计带通滤波器
Wn = [f_low, f_high] / (Fs/2);
if any(Wn >= 1)
Wn = min(0.99, Wn); % 防止归一化频率≥1
end
% 动态确定滤波器阶数
max_order = floor((segment_length-1)/3);
filter_order = min(max(1, max_order), 12); % 阶数限制在1-12
% 尝试filtfilt
try
[b, a] = butter(filter_order, Wn, 'bandpass');
filtered_noise = filtfilt(b, a, noise);
catch
try
b = fir1(filter_order*4, Wn, 'bandpass');
filtered_noise = filtfilt(b, 1, noise);
catch
filtered_noise = freq_domain_filter(noise, Fs, f_low, f_high);
end
end
% 应用渐变窗
fade_samples = min(50, floor(segment_length/10));
if fade_samples > 2
fade_in = linspace(0, 1, fade_samples);
fade_out = linspace(1, 0, fade_samples);
filtered_noise(1:fade_samples) = filtered_noise(1:fade_samples) .* fade_in;
filtered_noise(end-fade_samples+1:end) = filtered_noise(end-fade_samples+1:end) .* fade_out;
end
% 放入信号
signal(valid_idx) = filtered_noise;
% 归一化
if max(abs(signal)) > 0
signal = signal / max(abs(signal));
end
end
%% 辅助函数:生成单频信号
function [signal, t] = generate_single_tone(t_start, t_end, freq, Fs)
duration = t_end - t_start;
t = t_start:1/Fs:t_end;
signal = sin(2*pi*freq*t);
% 应用渐变窗
fade_samples = min(100, floor(length(signal)/4));
if fade_samples > 2
fade_in = linspace(0, 1, fade_samples);
fade_out = linspace(1, 0, fade_samples);
signal(1:fade_samples) = signal(1:fade_samples) .* fade_in;
signal(end-fade_samples+1:end) = signal(end-fade_samples+1:end) .* fade_out;
end
end
%% 辅助函数:生成脉冲信号
function [signal, t] = generate_impulse(t_center, f_low, f_high, Fs)
% 生成带限脉冲
center_freq = (f_low + f_high)/2;
bandwidth = f_high - f_low;
% 脉冲持续时间 (反比于带宽)
pulse_duration = min(0.1, 2/bandwidth);
t_start = t_center - pulse_duration/2;
t_end = t_center + pulse_duration/2;
t = t_start:1/Fs:t_end;
% 创建带限脉冲
pulse = sin(2*pi*center_freq*t) .* gausswin(length(t))';
% 归一化
signal = pulse / max(abs(pulse));
end
%% 辅助函数:频域滤波
function filtered = freq_domain_filter(signal, Fs, f_low, f_high)
N = length(signal);
freq = (0:N-1)*(Fs/N);
% 创建频域滤波器
H = ones(1, N);
low_idx = floor(f_low/(Fs/N)) + 1;
high_idx = ceil(f_high/(Fs/N)) + 1;
% 设置带通范围
H(1:low_idx) = 0;
H(high_idx:end) = 0;
H(N-high_idx+2:N) = 0; % 处理负频率部分
H(N-low_idx+2:N) = 0;
% 应用频域滤波
spectrum = fft(signal);
filtered_spectrum = spectrum .* H;
filtered = real(ifft(filtered_spectrum));
% 确保实信号
if isreal(signal)
filtered = real(filtered);
end
end
%% 修复后的可视化子函数 (完全兼容MATLAB 2019a)
function visualize_spectrogram(signal, Fs, grayImg, f1, f2, t1, t2)
fig = figure('Color', 'white', 'Position', [100, 100, 900, 700]);
try
% 时频分析参数
window = hamming(256);
noverlap = 200;
nfft = 1024;
% 计算时频谱
[s, freq, time] = spectrogram(signal, window, noverlap, nfft, Fs, 'yaxis');
% === 修复1: 正确获取waterfall对象句柄 ===
subplot(211)
h_waterfall = waterfall(time, freq, 10*log10(abs(s)));
% === 修复2: 设置曲面属性而非坐标轴属性 ===
set(h_waterfall, 'EdgeColor', 'interp', 'LineWidth', 0.5);
view([30, 60])
title('时频域瀑布图 (基于二值图像生成)')
xlabel('时间 (s)'), ylabel('频率 (Hz)'), zlabel('幅度 (dB)')
colormap jet
colorbar
grid on
ylim([f1, f2])
zlim([-50, 0])
% === 修复3: 兼容性图像显示 ===
subplot(212)
% 计算图像显示位置
time_vec = linspace(t1, t2, size(grayImg, 2));
freq_vec = linspace(f1, f2, size(grayImg, 1));
% 使用imagesc代替image确保兼容性
imagesc(time_vec, freq_vec, grayImg);
set(gca, 'YDir', 'normal')
title('输入的二值时频结构')
xlabel('时间 (s)'), ylabel('频率 (Hz)')
colormap(gca, gray(256)) % 明确指定灰度图
colorbar
axis tight
sgtitle(sprintf('时频结构可视化: %.1f-%.1fHz, %.1f-%.1fs', f1, f2, t1, t2), 'FontSize', 14)
catch ME
warning('可视化失败: %s', ME.message);
% 备用绘图方案
subplot(211)
plot((0:length(signal)-1)/Fs, signal);
title('生成的时域信号');
xlabel('时间 (s)'), ylabel('幅度');
grid on;
subplot(212)
imagesc(grayImg);
title('输入的二值时频结构');
colormap gray;
colorbar;
end
end
请根据用户需求修改生成新的函数,generate_tf_waterfall中引入一个新的变量n,其功能如下:
变量n表示将信号在频域上切分为均等的n段,每一段由一个独立的信号源实现切分后的瀑布图效果,返回值composite_signal修改为一个2维数组,表示每个信号源输出的时域信号,visualize_spectrogram中绘制的瀑布图应为各信号源信号合成后的瀑布图
最新发布