活动介绍
file-type

使用MATLAB绘制10Hz正弦信号频谱图教程

版权申诉

RAR文件

2KB | 更新于2024-12-07 | 189 浏览量 | 0 下载量 举报 收藏
download 限时特惠:#14.90
这个过程涉及到信号处理的基本概念,以及Matlab在信号分析中的具体应用方法。" 知识点详细说明: 1. 正弦信号的概念: 正弦信号是物理学和信号处理中的一种基本信号,其表达式通常为 y(t) = A * sin(ωt + φ),其中 A 代表振幅,ω 代表角频率,φ 代表相位。角频率 ω 与频率 f 的关系为 ω = 2πf。因此,一个频率为10Hz的正弦信号意味着其周期 T = 1/f = 0.1秒。 2. 功率定义: 在信号处理中,信号的功率通常表示为信号强度的平方的均值,对于正弦信号来说,功率 P 可以表示为 P = A^2/2。本例中给出功率为1,意味着振幅 A 的平方为2,因此振幅 A = sqrt(2)。 3. Matlab编程基础: Matlab是一种用于算法开发、数据可视化、数据分析以及数值计算的高级编程语言和交互式环境。在本资源中,Matlab将被用来生成正弦信号并绘制其频谱。 4. 信号的采样与离散傅里叶变换(DFT): 在数字信号处理中,连续信号需要被采样以转换为离散信号。采样定理指出,采样频率必须至少是信号最高频率的两倍,以避免混叠现象。对于10Hz的信号,如果采样率没有给出,则一般默认使用Matlab的默认采样频率(通常为1000Hz以上)。Matlab使用离散傅里叶变换(DFT)算法计算信号的频谱,该算法将时域信号转换为频域信号。 5. Matlab中的频谱绘制: 频谱图是信号频域表示的一种直观展示,显示了信号在不同频率下的强度。在Matlab中,可以使用快速傅里叶变换(FFT)函数来高效计算DFT,然后使用绘图函数如plot来绘制频谱图。Matlab内置函数fft可以用来计算信号的FFT。 6. 使用Matlab编程实现: 首先,在Matlab中使用sin函数生成10Hz的正弦波信号。然后,通过适当采样率对信号进行采样,以保证信号的正确表示。接着,使用fft函数计算该离散时间信号的频谱,并通过plot函数绘制其幅度谱。 7. 分析和解释绘制结果: 绘制出的频谱图会显示出一个峰值,位于10Hz处,表示信号的主要频率成分。频谱图上其他位置的数值大小可以表示噪声或其他频率成分的强度。如果信号中只包含10Hz的单一频率,频谱图中应该只看到一个明显的尖峰。 8. 文档mat sin.doc的使用: 虽然提供的文件列表中只有一个名为"mat sin.doc"的文件,但是我们可以合理推测该文档应包含上述所有概念的详细解释和Matlab代码实现步骤,可能是用于指导用户完成整个过程的教程文档。 总结上述知识点,本资源主要围绕在Matlab环境下创建特定频率和功率的正弦信号,并使用Matlab内置工具绘制该信号的频谱图。这不仅涉及到了信号处理的基本理论,还包括了Matlab在工程实践中的应用,对于理解和掌握数字信号处理的基本方法具有重要的实践意义。

相关推荐

filetype

现在已有两个函数如下: %% 修改后的 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中绘制的瀑布图应为各信号源信号合成后的瀑布图

filetype

%% 1. 生成合成甲烷信号和多个含噪信号 close all; % 关闭所有图形窗口 clear all; % 清除所有变量 clc; % 清除命令窗口 Fs = 1000; % 采样频率 t = 0:1/Fs:0.999; % 时间向量 clean_signal = sin(2*pi*5*t); % 清洁信号:5 Hz正弦波 snr_dB = -5; % 目标信噪比,单位dB N = 5; % 训练含噪信号的数量 training_noisy_signals = cell(1,N); % 初始化训练含噪信号的单元数组 for i=1:N training_noisy_signals{i} = awgn(clean_signal, snr_dB, 'measured'); % 生成含噪信号 end % 生成一个测试含噪信号 test_noisy_signal = awgn(clean_signal, snr_dB, 'measured'); % 生成测试含噪信号 %% 2. 数据预处理(无重叠分帧) window_size = 40; % 帧大小 overlap = 0; % 无重叠 step = window_size - overlap; % 步长 % 对清洁信号分帧 T = buffer(clean_signal, window_size, overlap, 'nodelay'); % 分帧清洁信号 % 对训练含噪信号分帧并拼接 X_train = []; % 初始化训练输入矩阵 for i=1:N X_i = buffer(training_noisy_signals{i}, window_size, overlap, 'nodelay'); % 分帧每个训练含噪信号 X_train = [X_train, X_i]; % 水平拼接 end T_train = repmat(T, 1, N); % 重复T N次 % 随机打乱训练数据以提高泛化能力 num_frames_total = size(X_train, 2); % 总帧数 indices = randperm(num_frames_total); % 随机排列索引 X_train = X_train(:, indices); % 打乱训练输入 T_train = T_train(:, indices); % 打乱训练目标 %% 3. 创建并训练BP神经网络 net = feedforwardnet([40 20], 'trainbr'); % 两隐藏层:40和20个神经元 net.divideFcn = 'divideblock'; % 块划分 net.divideParam.trainRatio = 70/100; % 70%训练 net.divideParam.valRatio = 15/100; % 15%验证 net.divideParam.testRatio = 15/100; % 15%测试 net.performParam.regularization = 0.01; % L2正则化 net.trainParam.epochs = 1000; % 最大迭代次数 net.trainParam.lr = 0.02; % 学习率 net.trainParam.showWindow = true; % 显示训练窗口 % 训练网络 [net, tr] = train(net, X_train, T_train); % 训练网络 save('denoise_net_v2.mat', 'net'); % 保存训练好的网络 %% 4. 测试降噪效果 load('denoise_net_v2.mat'); % 加载训练好的网络 % 对测试含噪信号分帧 X_test = buffer(test_noisy_signal, window_size, overlap, 'nodelay'); % 分帧测试含噪信号 % 使用网络进行降噪 denoised_signal = sim(net, X_test); % 降噪帧 % 重构完整的降噪信号 denoised_full_test = denoised_signal(:)'; % 转换为行向量 %% 5. 可视化对比 figure(1); plot(t, clean_signal, 'b', 'LineWidth', 1.5); hold on; % 绘制清洁信号,蓝色 plot(t, test_noisy_signal, 'g', 'LineWidth', 1); % 绘制测试含噪信号,绿色 plot(t, denoised_full_test, 'r', 'LineWidth', 1.5); % 绘制降噪信号,红色 xlabel('Time (s)'); ylabel('Amplitude'); % 标注轴 legend('Clean Signal', 'Noisy Signal', 'Denoised Signal'); % 添加图例 title('BP Neural Network Denoising Effect Comparison'); % 添加标题 figure(2); plot(t, denoised_full_test, 'r', 'LineWidth', 2); % 绘制降噪信号 xlabel('Time (s)'); ylabel('Amplitude'); % 标注轴 %% 6. 性能评估 original_snr = 10*log10(var(clean_signal)/var(test_noisy_signal - clean_signal)); % 计算测试信号的原始SNR denoised_snr = 10*log10(var(clean_signal)/var(denoised_full_test - clean_signal)); % 计算降噪后的SNR fprintf('Original SNR: %.2f dB\n', original_snr); % 打印原始SNR fprintf('Denoised SNR: %.2f dB\n', denoised_snr); % 打印降噪后的SNR fprintf('SNR improvement: %.2f dB\n', denoised_snr - original_snr); % 打印SNR提升量这是目前采用的BP神经网络信号降噪处理算法