赞
踩
clear; clc; close all;
% 生成一个序列,用希尔伯特分析信号
xr = [1 2 3 4];
x = hilbert(xr)
% 其虚部
x_im = imag(x)
% 其实部
x_re = real(x)
% 快速傅里叶变换
x_dft = fft(x)
xr_dft = fft(xr)
x = 1.0000 + 1.0000i 2.0000 - 1.0000i 3.0000 - 1.0000i 4.0000 + 1.0000i x_im = 1 -1 -1 1 x_re = 1 2 3 4 x_dft = 10.0000 + 0.0000i -4.0000 + 4.0000i -2.0000 + 0.0000i 0.0000 + 0.0000i xr_dft = 10.0000 + 0.0000i -2.0000 + 2.0000i -2.0000 + 0.0000i -2.0000 - 2.0000i
% 希尔伯特函数为有限长数据找到精确的解析信号 % 也可以通过使用FIR希尔伯特变换滤波器来生成解析信号,以计算虚部的近似 % 生成一组频率为203,721,1001Hz的正弦序列 % 采样频率为10kHz,时长为1秒 clear; clc; close all; fs = 1e4; t = 0:1/fs:1; x = 2.5 + cos(2*pi*203*t) + sin(2*pi*721*t) + cos(2*pi*1001*t); y = hilbert(x); plot(t,real(y),t,imag(y)) xlim([0.01 0.03]) legend('real','imaginary') title('hilbert Function') xlabel('Time (s)') % 计算原始序列和解析信号的Welch估计功率谱密度 % 窗长为256,无重叠部分的汉明窗 % 验证解析信号在负频率处没有功率 pwelch([x;y].',256,0,[],fs,'centered') % help pwelch legend('Original','hilbert') % 设计一个60阶希尔伯特变换FIR滤波器,过渡带宽为400Hz fo = 60; d = designfilt('hilbertfir','FilterOrder',fo, ... 'TransitionWidth',400,'SampleRate',fs); freqz(d,1024,fs) % 对正弦序列进行滤波以逼近解析信号的虚部 hb = filter(d,x) % FIR Filter grd = fo/2; y2 = x(1:end-grd) + 1j*hb(grd+1:end); t2 = t(1:end-grd); plot(t2,real(y2),t2,imag(y2)) xlim([0.01 0.03]) legend('real','imaginary') title('FIR Filter') xlabel('Time (s)') % 估计近似解析信号的功率谱密度(PSD),并与希尔伯特结果进行比较 pwelch([y;[y2 zeros(1,grd)]].',256,0,[],fs,'centered') legend('hilbert','FIR Filter')
上图为解析信号x(t)的实部与虚部;
计算原始序列和解析信号的Welch估计功率谱密度;解析信号频率大于0
设计一个60阶希尔伯特变换FIR滤波器,过渡带宽为400Hz;
对比前面的图,对正弦序列进行滤波也可以逼近解析信号的虚部;
估计近似解析信号的功率谱密度(PSD),并与希尔伯特结果进行比较,与前图效果相差不大。
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。