希尔伯特
最重要的就两句
z=hilbert(x) am2=abs(z)
首先,不是第一句完成了希尔伯特变换就是包络了,它只是完成了希尔伯特变换,相位偏移了90度,而且要使用abs()而不是imag(),下面代码运行就可以得到对应的图,但是前提是要把数据包添加到matlab的文件夹里
链接:http://pan.baidu.com/s/1kU9DuAF 密码:gblw
,并且下列代码的m文件要和数据包在同一路径下
clc;
clear;
patterns = [];
targets = [];
Fs=50e6;%采样率
%Abstract feature
for i=1:1
floder = ['IQ_sincos_f2.395000e+09_s5.000000e+07_Tx' num2str(i) '\data'];
fea_train =[];
for j=1:1
path = [floder '\IQ_sincos_f2.395000e+09_s5.000000e+07_Tx' num2str(i) '_' num2str(j) '.pcm'];
fileId = fopen(path,'r');
x = fread(fileId,500,'float32');
lengthofsignal=length(x);
T=1/Fs;
t=(0:lengthofsignal-1)*T; %对应的时间序列
%绘制原始图像
figure();
subplot(211)
plot(t,x,'b');
grid on
xlabel('时间');
ylabel('幅度')
hold on
z=hilbert(x);
am=abs(z);
plot(t,am,'r')
title('希尔伯特变换后使用abs函数');
legend('蓝色是原始波形','红色是abs求出的包络')
%提取信号包络并画图
subplot(212)
plot(t,x,'b');
grid on
xlabel('时间');
ylabel('幅度')
hold on
%提取信号包络并画图,法二
z=hilbert(x);
yi = imag(z);
plot(t,yi,'g')
title('希尔伯特变换后使用imag函数');
legend('绿色是imag求出的')
end
end
其次,使用希尔伯特可以得到包络,但是得到包络的方法不止一只,还可以使用低通滤波,主程序如下:
clc;
clear;
patterns = [];
targets = [];
Fs=12000;%采样率
%Abstract feature
for i=1:1 %只用了一个来测试,如果要用多个改掉后面的1
floder = ['IQ_sincos_f2.395000e+09_s5.000000e+07_Tx' num2str(i) '\data'];
fea_train =[];
for j=1:1
path = [floder '\IQ_sincos_f2.395000e+09_s5.000000e+07_Tx' num2str(i) '_' num2str(j) '.pcm'];
fileId = fopen(path,'r');
x = fread(fileId,10240,'float32');
lengthofsignal=length(x);
T=1/Fs;
t=(0:lengthofsignal-1)*T; %对应的时间序列
%绘制原始图像
figure(1);
plot(t,x);
grid on
title(['TX',num2str(i),'原始图像并采用希尔伯特求包络']);
xlabel('Time');
ylabel('Amp')
hold on
z=hilbert(x);%做希尔伯特变换
%am1=sqrt(abs(x).^2+abs(z).^2);%使用此种函数处理后的包络没有用abs的好
am2=abs(z);
%plot(t,am1,'r')
plot(t,am2,'g')
axis([0 0.25 -1 1.5]); % 设置坐标轴在指定的区间
figure(2)
plot(t,x);
hold on
re=envelope(x,Fs)
plot(t,re,'r--');title('低通滤波求包络');
axis([0 0.25 -1 1.5]); % 设置坐标轴在指定的区间
figure(3)
plot(t,am2,'g')
hold on
plot(t,re,'r--')
title('对比')
legend('绿色是希尔伯特求出的包络','红色是低通滤波求出的包络')%文本标注
end
end
其中调用的envelope.m文件
function y=envelope(signal,Fs)
%Envelope Detection based on Low pass filter and then FFT
[a,b]=butter(2,0.1);%butterworth Filter of 2 poles and Wn=0.1
%sig_abs=abs(signal); % Can be used instead of squaring, then filtering and
%then taking square root
sig_sq=2*signal.*signal;% squaring for rectifing
%gain of 2 for maintianing the same energy in the output
y_sq = filter(a,b,sig_sq); %applying LPF
y=sqrt(y_sq);%taking Square root
%advantages of taking square and then Square root rather than abs, brings
%out some hidden information more efficiently
N=3000;T=N/Fs;
sig_f=abs(fft(y(1:N)',N));
sig_n=sig_f/(norm(sig_f));
freq_s=(0:N-1)/T;
从图上可以看出,采用低通滤波的包络更为集中,更稳定一些