首页 > 分享 > 小波变换与信号处理：去噪方法和频谱分析

小波变换与信号处理：去噪方法和频谱分析

花匠小妙招
2025-02-08 10:17

傅里叶变换

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

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

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

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

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

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

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

ψ ∗ ( a , b ) = 1 a ψ ( t − b a ) psi ^*(a,b) = frac{1}{sqrt{a}} psi(frac{t-b}{a}) ψ∗(a,b)=a

1ψ(at−b)

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

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

在这里插入图片描述

分解层数的选择

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

例如：
信号持续时间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');

12345678910111213141516171819202122

二者相减即为水面振动信号
在这里插入图片描述

小波去噪

为什么要使用阈值

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

阈值获取

阈值获取的函数有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'硬阈值，大于阈值的点保持不变。一般来说，硬阈值处理后的信号比用软阈值处理后的信号更粗糙。 1234

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

层数选择

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

在这里插入图片描述

上图是指厘米量级海面干扰下小波分解得到几十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('周期图法频率估计');

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125