一维信号去噪分析

一维信号滤波分析

读取实验数据，绘制一维时域谱，如图1所示

23.8

23.8523.923.95

2424.0524.124.1524.2

0.0408

0.0410.04120.04140.04160.0418

0.0420.0422

0.0424运动位置

信号幅值

原始信号

图1

下面对信号进行FFT 变换： Fs=1024;%假定采样频率 N=597;%实际采样点数

y=xlsread('C:\Documents and Settings\Administrator\桌面

\Book2.xls'); y=y';

y1=fft(y);

y2=10*log10((abs(y1).^2)); figure; plot(y2);

title('fft 变换')%做FFT 变换 Ay2=abs(y2); %取模

Ay2=Ay2/(N/2); %换算成实际的幅度 Ay2(1)=Ay2(1)/2;

F=([1:N]-1)*Fs/N; %换算成实际的频率值，Fn=(n-1)*Fs/N figure

plot(F(1:N/2),Ay2(1:N/2)); %显示换算后的FFT 模值结果 title('幅度-频率曲线图');

得到信号频谱如图2

0100200300400500600

0.05

0.1

0.15

0.2

0.25

0.3

0.35

幅度-频率曲线图

图2

由于大多数信号噪声频率都为高频，所以可以先尝试设计一个简单的带通滤波器，这里采用窗函数法，根据FFT 变换得到的频谱图取通带为50~350Hz ，算法函数如下： Fs=1024; N=597; n=36;

y=xlsread('C:\Documents and

Settings\Administrator\

桌

面

\Book2.xls'); y=y';

x=xlsread('C:\Documents and Settings\Administrator\桌面

\Book1.xls'); x=x';

Wn=[50,350]*2/Fs;%归一化频率 b=fir1(n,Wn);

yf=filter(b,1,y); figure; plot(x,y);

xlabel('运动位置'); ylabel('信号幅值');

title('窗函数法滤波后的信号')

滤波器频谱响应如图3

00.10.2

0.30.40.50.60.70.80.91

-3000

-2000

-1000

Normalized Frequency (?π rad/sample)

P h a s e (d e g r e e s )

00.10.2

0.30.40.50.60.70.80.91

-100

-50

50

Normalized Frequency (?π rad/sample)

M a g n i t u d e (d B )

滤波后得到滤波信号如图4

23.8

23.8523.923.95

2424.0524.124.1524.2

0.0408

0.0410.04120.04140.04160.0418

0.0420.0422

0.0424运动位置

信号幅值

窗函数法滤波后的信号

图4

从图4可以看出，噪声部分被滤除，但这样设计的带通滤波器对去噪没有明显的效果，并且在未知采样频率的情况下对不确定信号进行FFT 变换的效果并不理想，而小波分析则克服了这种缺陷，它对随机信号具有多分辨率分析的特点，在时域和频域都有表征信号局部信息的能力，时间窗和频率窗都可以根据信号的具体形态动态调整，在一般情况下，在低频

部分可以采用较低的时间分辨率，而提高频率的分辨率，在高频情况下可以用较低的频率分辨率来换取精确的时间定位。 针对以上实验数据，这里采用的思路是先对信号进行低层分解强度去噪，然后通过阈值信号处理，得到最终的滤波信号，小波算法函数设计如下： y=xlsread('C:\Documents and Settings\Administrator\桌面\Book2.xls'); y=y';

x=xlsread('C:\Documents and

Settings\Administrator\

桌

面

\Book1.xls'); x=x';

[c,l]=wavedec(y,1,'db1');%用db1小波对原始信号进行1层分解并提取系数 ca1=appcoef(c,l,'db1',1);%提取近似系数 cd1=detcoef(c,l,1);%提取1阶详细系数

cdd1=zeros(1,length(cd1));%对信号进行强制性去噪处理并显示结果 c1=[ca1 cdd1];

s1=waverec(c1,l,'db1');%多层次的一维小波重构

N=2;%N 越大，分解的层数越多，则能滤除更多噪声，同时取值过大也会导致信号失真 yd=wden(s1,'rigrsure','h','mln',N,'sym20');%rigrsure 阈值信号处理 figure;

plot(x,yd);

xlabel('运动位置'); ylabel('信号幅值');

title('小波分析滤波后的信号')

得到滤波信号如图5

23.8

23.8523.923.95

2424.0524.124.1524.2

0.0408

0.0410.04120.04140.04160.0418

0.0420.0422

0.0424运动位置

信号幅值

小波分析后滤波信号

图5

如图5所示，信号有明显改善，大部分噪声已去除，但是波形还是很不规则。假设信号是

由多个单频信号组成的组合信号，我们可以针对信号的频谱直接设计一个滤波器，而通过对信号进行FFT 变换并不能得到想要的结果，所以这里采用另外一种谱分析方法——功率谱估计，它是专门针对随机信号所采取的谱分析方法。功率谱估计分为非参数方法、参数方法、子空间方法三大类，子空间方法相对来说有更高的分辨率，它又可分为MUSIC 方法和特征向量法。其中，特征向量法主要适用于混有高斯白噪声的正弦信号的功率谱估计，而MUSIC 法适合更为普遍情况下正弦信号的功率谱估计，所以这里采用MUSIC 法对信号进行初步分析，得到功率谱图如图6

00.10.2

0.30.40.50.60.70.80.91

-50

50

100

150

200

250

X: 0.1641Y: 79.83

Normalized Frequency ( rad/sample)

P o w e r (d B )

Pseudospectrum Estimate v ia MUSIC

图6

由图6可以看出，信号在频率为0.16，0.62两点处分别出现峰值，我们假设信号就是由这两种频率构成的单频信号的组合，针对此信号频率设计一个滤波器，其基本步骤如下： n=50;

f=[0 0.12 0.13 0.20 0.21 0.58 0.59 0.65 0.66 1]; m=[0 0 1 1 0 0 1 1 0 0]; [b a]=firls(n,f,m); 滤波器频率响应如图7

0.1

0.2

0.30.40.50.60.70.80.9

1

-1500

-1000-5000

500

Normalized Frequency (?π rad/sample)

P h a s e (d e g r e e s )

00.10.2

0.30.40.50.60.70.80.91

-100

-50

50

Normalized Frequency (?π rad/sample)

M a g n i t u d e (d B )

图7 先对用该滤波器对信号进行初步滤波，然后用平滑函数进行处理，最后用前面的小波分析得到滤波信号，其算法如下：

y=xlsread('C:\Documents and Settings\Administrator\桌面\Book2.xls'); y=y';

x=xlsread('C:\Documents and Settings\Administrator\桌

面

\Book1.xls'); x=x' n=90;

f=[0 0.12 0.13 0.20 0.21 0.58 0.59 0.65 0.66 1]; m=[0 0 1 1 0 0 1 1 0 0]; [b a]=firls(n,f,m); y1=filter(b,a,y);

y1(1,[1:100])=y1(1,[371:470]); y2=smooth(y1);

y3=wden(y2,'rigrsure','h','mln',2,'sym20'); figure;

plot(x,y3); xlabel('运动位置'); ylabel('信号幅值'); title('滤波信号')

得到的滤波信号如图8

23.8

23.85

23.9

23.95

2424.05

24.1

24.15

24.2

0.005

0.010.0150.020.0250.030.0350.040.045运动位置

信号幅值

滤波信号

图8 由此得到的的滤波信号虽然规则平滑，可是信号幅值与实验数据相比有较大差距，所以此种方法的可行性有待进一步研究

此外，也可以采用LMS 算法进行自适应滤波，得到的滤波信号也比较理想，其算法步骤如下：

num=597;%采样点数 delay=10;%延迟间隔 N=10;%滤波器抽头数 u=0.001;%调整步长 k=1:num;

y=xlsread('C:\Documents and

Settings\Administrator\

桌

面

\Book2.xls'); yf=y';

x=xlsread('C:\Documents and Settings\Administrator\桌面

\Book1.xls'); x=x';

yf_in=[yf zeros(1,delay)];%当前的输入信号 yf_delay=[zeros(1,delay),yf];%延迟的输入信号 M=num+delay;

y2=zeros(1,M);%输出初始值 w=zeros(1,N);%滤波系数初始值 for n=N:M-N+1

yfyf=yf_delay(n+N-1:-1:n);%滤波器输入值 y2(n)=w*yfyf';%滤波器输出值 e(n)=yf_in(n);

w=w+2*u.*e(n).*yfyf;%系数调整 end

figure;

y3=e(10:num);%提取信号

y4=wden(y3,'rigrsure','h','mln',2,'sym20');%小波进一步滤波 y4=smooth(y4);

x=x(:,[1:588]);%保持向量维数一致 plot(x,y4);

xlabel('运动位置'); ylabel('信号幅值');

title('LMS 算法滤波后的信号');

得到的处理结果如图9

23.8

23.8523.923.95

2424.0524.124.1524.2

0.041

0.0412

0.0414

0.0416

0.0418

0.042

0.0422

运动位置

信号幅值

LMS 算法滤波后的信号

图9

clear all

close all

clc

V1 = rand(20,1)*700+200;

subplot 231; plot(V1); ylim([200 700]); grid;

y=mean5_3(V1,5);

subplot 232; plot(y); ylim([200 700]); grid;

y=mean5_3(V1,10);

subplot 233; plot(y); ylim([000 700]); grid;

y=mean5_3(V1,50);

subplot 234; plot(y); ylim([200 700]); grid;

y=mean5_3(V1,200);

subplot 235; plot(y); ylim([200 700]); grid;

y=mean5_3(V1,500);

subplot 236; plot(y); ylim([200 700]); grid;

%------------------------------------------------------------ function y = mean5_3(x, m)

% x为被处理的数据

% m 为循环次数

n=length(x);

a=x;

for k=1: m

b(1) = (69*a(1) +4*(a(2) +a(4)) -6*a(3) -a(5)) /70;

b(2) = (2* (a(1) +a(5)) +27*a(2) +12*a(3) -8*a(4)) /35;

for j=3:n-2

b (j) = (-3*(a(j-2) +a(j+2)) +12*(a(j-1) +a(j+1)) +17*a(j)) /35;

end

b (n-1) = (2*(a(n) +a(n-4)) +27*a(n-1) +12*a(n-2) -8*a(n-3)) /35;

b (n) = (69*a(n) +4* (a(n-1) +a(n-3)) -6*a(n-2) -a(n-4)) /70;

a=b;

end

y =a;

subplot 231;

plot(V1);

grid;

y=mean5_3(V1,5);

subplot 232; plot(y);

grid;

y=mean5_3(V1,10);

subplot 233;

plot(y);

grid;

y=mean5_3(V1,50);

subplot 234; plot(y);

grid;

y=mean5_3(V1,200);

subplot 235;

plot(y);

grid;

y=mean5_3(V1,500);

subplot 236; plot(y);

grid;