小波变换与小波阈值去噪

傅里叶变换

三角函数基，缩得窄对应高频，伸得宽对应低频；基函数不断和信号相乘，某个尺度乘出来的结果，就可以理解成信号所包含的当前尺度对应频率成分有多少。于是，基函数会在某些尺度下，与信号相乘得到一个很大的值，因为此时二者有一种重合关系。那么我们就知道信号包含该频率的成分的多少。

傅里叶变换可以分析信号的频谱
但对于非平稳过程（频率随时间变化），傅里叶变换有局限性
他只能获取一段信号总体上包含的频率成分，但无法得到个频率成分出现的时刻。
因此是与相差很大的信号，可能频谱图一样。

短时傅里叶变换

把整个时域过程分成无数个等长的小过程，每个小过程近似平稳，再傅里叶变换，就知道在哪个时间点上出现什么频率。但STFT仍存在缺陷，窗太窄，窗内信号太短，导致频率分袖析够精准，频率分辨率差；窗太宽，时间分辨率差。

对于时变的非稳态信号，高频适合小窗口，低频适合大窗口，然而STFT的窗口是固定的。

小波变换

小波变换和STFT的出发点不同，STFT是给信号加窗，分段做FFT；小波变换将无限长的三角函数基换成$\mathbf{有限长会衰减的小波基}$。

 小波母函数满足紧支撑性、波动性、容许条件、正交性。

（1）紧支撑性：仅在一小部分定义域里不为零。具有紧支撑性的基函数，在原信号的时间轴上平移，就相当于对原信号加窗操作。
（2）波动性：在定义域内积分值为0；
（3）容许条件和正交性使变换可逆。

$${\psi }^{\ast }(a,b)=\frac{1}{\sqrt{a}}\psi (\frac{t-b}{a})$$

a较小，相当于挤压，提高了频率，窗子变小；a较大，相当于拉伸，降低了频率，窗子变大。
低频，宽窗，好的频域分辨率；高频，窄窗，好的时间分辨率。

小波阈值法进行信号去噪，对信号进行五层小波分解，每进行一层分解时都将该层得到的细节系数进行一次阈值化处理。

分解层数的选择

分解层数越大，噪声和信号表现得不同特性越明显，越有利于二者的分离；但重构到的信号失真也比较大，一定程度上又会影响最终去噪效果。

例如：
信号持续时间6秒，30000个点，采样率为5000Hz，由采样定理知，该信号的最大频率为2500Hz，对该信号做3层DWT，一阶细节的频段为1250-2500Hz，一阶近似的频段小于1250Hz，二阶近似的频段小于625Hz，三阶近似的频段小于312.5Hz，四阶近似的频段小于156.25Hz，五阶近似的频段小于78.125Hz，

[c,l]=wavedec(y_inf,5,'db4'); 
%取第5层低频近似系数
ca5=appcoef(c,l,'db4',5);
%取各层高频细节系数
cd5=detcoef(c,l,5);
cd4=detcoef(c,l,4);
cd3=detcoef(c,l,3);
cd2=detcoef(c,l,2);
cd1=detcoef(c,l,1);
thr=thselect(y_inf,'sqtwolog'); % 阈值获取
%进行软阈值处理
ysoft5=wthresh(cd5,'s',thr);
ysoft4=wthresh(cd4,'s',thr);
ysoft3=wthresh(cd3,'s',thr);
ysoft2=wthresh(cd2,'s',thr);
ysoft1=wthresh(cd1,'s',thr);
c1=[ca5;ysoft5;ysoft4;ysoft3;ysoft2;ysoft1];
Y=waverec(c1,l,'db4');
y_inf_2 = y_inf - Y;
figure;
subplot(1,2,1);plot(tau,y_inf);xlabel('t/s');ylabel('幅值');title('含噪信号y_inf');
subplot(1,2,2);plot(tau,Y);xlabel('t/s');ylabel('幅值');title('去噪信号Y');

二者相减即为水面振动信号

小波去噪

为什么要使用阈值

通常从设备上采集到的信号是有效信号和噪声的和。
在小波域，有效信号的对应的系数很大，而噪声对应的系数很小。

阈值获取

阈值获取的函数有ddencmp、thselect 、wbmpen 、wwdcbm

thr = thselect(X,TPTR)
TPTR=‘rigrsure’ ‘sqtwolog’ ‘minimaxi’

阈值去噪

实现阈值去噪的函数有wden wthresh wthcoef

y = wthresh(X, SORH,thr);
SORH = 's'  'h'
's'软阈值，信号绝对值与阈值进行比较，小于或等于阈值的点变为零，大于阈值的点为该点与阈值的差值；
'h'硬阈值，大于阈值的点保持不变。一般来说，硬阈值处理后的信号比用软阈值处理后的信号更粗糙。

用软阈值函数是为了解决硬阈值函数“一刀切”导致的影响。

层数选择

小波变换只对信号的低频部分作进一步分解，对高频不再继续分解。所以小波变换能够很小的表征以低频信息为主要成分的信号。
小波层数越多会丢失细节，层次越少不能去除噪声。
高层分解的小波系数对应的是低频部分。

上图是指厘米量级海面干扰下小波分解得到几十nm量级的水面微幅振动。
原理是：海面干扰低频，信号/声源高频，小波去噪重构得到的低频信号，然后相减即为信号。

%太赫兹雷达，发射chirp，接收，提取相位

close all;clear
f = 200; %水下声源频率
P = 170;    %点声压
omega = 2*pi*f; %角频率:/s
rho = 1e3;
c = 1500;   %声波在水介质的传播速度
lamda = c/f ;
k = 2*pi/lamda; %水面波动的波数:/m
mu = 1;%水介质粘度：Pa s
g = 9.8;%重力加速度:N/kg
sigma = 72.75e-3;%表面张力系数:N/m
alpha = 4*mu*k^2*sqrt(k*g+sigma*k^3/rho)/(rho*g+3*sigma*k^2/rho);
ymax = 2*P/omega/rho/c;
%-----------------------------雷达---------------
h = 2;%雷达水上高度
fc = 122.5e9;%雷达中心频率,初始频率为122GHz

c0 =  3e8; %光速
lamda_r = c0/fc;
Br = 1e9;%雷达带宽
dr = c0/2/Br;   %距离分辨率
disp('距离分辨率（cm）：')
disp(dr*100 )
Tp= 0.1e-3; %chirp持续时间
Tr = 0.2e-3; %扫频周期
Kr = Br/Tp;
fs = 1e6;
pulse_len = round(Tp*fs);   %脉冲长度
xx = 0.5; %声源处声致水面波振动
T = 0.6; %总时长6s
t = (0:pulse_len-1)/fs;
r = t*c0/2/1e3;
pulse_num =round(T/Tr) ;   %脉冲数
tau = (0:pulse_num-1)*Tr;
Mix_IF_fft = zeros(pulse_num, pulse_len);
%------------------加入厘米两级干扰
ht = 0.003*cos(-2*pi*100*linspace(0,0.1,pulse_num));
% ht = zeros(1,pulse_num);
% figure;plot(ht);title('海面干扰');
%----------------模拟雷达测距------------------------------
for i = 1:pulse_num
    ti = (i-1)*Tr + t;     %该脉冲时间

    Tx = exp(1j*2*pi*fc *t  + 1j*pi* Kr* t.^2);   %脉冲时间都是从零算的
    yt =  ymax * exp(-alpha * xx) * cos( - omega*ti);
    

    td = 2*(h-ht(i)-yt)/c0;
    Rx = exp(1j *2*pi*fc *(t-td) + 1j*pi* Kr*(t-td).^2 )  ;   %接收到的回波
    Mix = Tx .* conj(Rx);
    Mix_IF_fft(i,:)= fft(Mix, pulse_len);    %时间-距离谱图
end
% figure; imagesc(r,tau,abs(Mix_IF_fft));ylabel('时间/s');xlabel('距离/m');title('距离时间谱图');
%----------------选择正确的距离门----------------
tmp = zeros(1,pulse_len);
for i = 1:pulse_len
    tmp(1,i) = sum(abs(Mix_IF_fft(:,i)).^2)/pulse_num;
end
[~,indx] = max(tmp);
%---------------提取相位----------------------
phase = angle(Mix_IF_fft(:,indx) .* exp(-1j*2*pi*fc*2*h/c0));
y_inf = unwrap(phase)*lamda_r/4/pi;
figure;plot(tau,y_inf);title('推断出的水面振动幅度');
% hold on;plot(tau,hn);legend('推断出的水面振动幅度','加入的海面波动');hold off;

% --------------小波分解--------------------
[c,l]=wavedec(y_inf,4,'db4');
%取第5层低频近似系数
ca6=appcoef(c,l,'db4',4);
%取各层高频细节系数
% cd6=detcoef(c,l,6);
% cd7=detcoef(c,l,7);
% cd6=detcoef(c,l,6);
% cd5=detcoef(c,l,5);
cd4=detcoef(c,l,4);
cd3=detcoef(c,l,3);
cd2=detcoef(c,l,2);
cd1=detcoef(c,l,1);
% [thr,~] = ddencmp('den','wv',y_inf);
thr=thselect(y_inf,'sqtwolog'); % 阈值获取
%进行软阈值处理
% ysoft6=wthresh(cd6,'s',thr);
% ysoft7=wthresh(cd7,'s',thr);
% ysoft6=wthresh(cd6,'s',thr);
% ysoft5=wthresh(cd5,'s',thr);
ysoft4=wthresh(cd4,'s',thr);
ysoft3=wthresh(cd3,'s',thr);
ysoft2=wthresh(cd2,'s',thr);
ysoft1=wthresh(cd1,'s',thr);
c1=[ca6;ysoft4;ysoft3;ysoft2;ysoft1];
Y=waverec(c1,l,'db4');
%不断去除细节，重构得到低频分量Y
y_inf_2 =y_inf- Y;
figure;plot(tau,y_inf);xlabel('t/s');ylabel('幅值');title('含噪信号y_inf');
hold on;plot(tau,y_inf_2);hold off;

% ---------卡尔曼滤波--------
length_win= 30000;
num = length(y_inf_2)/length_win;
Xn = y_inf_2;
for k=1:num
    Z=y_inf_2((k-1)*length_win+1:k*length_win);
    r=5;
    X=Z(1);
    T=1;
    V=1;
    Q=T^2*V^2 ;
    R=r^2;
    P=X^2;
    A=1;
    H=1;
    tmp=KaermanFilter(X,Z, A, Q, H, R,P);
    Xn((k-1)*length_win+1:k*length_win) = tmp;
end
figure;plot(tau,y_inf_2);xlabel('t/s');ylabel('幅值');title('去噪信号');
hold on;plot(tau,Xn);hold off;
%--------------频率估计------------------
nfft = pulse_num;
Fs = 1/Tr;
% y_inf_2_hilt = hilbert(y_inf_2);
[P,f] = periodogram(hilbert(Xn),[],nfft,Fs);
% figure;plot(linspace(Fs/2,-Fs/2,length(y_inf)),P);title('周期图法频率估计');
figure;plot(f,P);title('周期图法频率估计');