一、简介
基于SVD(奇异值分解)的去噪声技术属于子空间算法的一种。简单的来说我们希望将带噪信号向量空间分解为分别由纯净信号主导和噪声信号主导的两个子空间,然后通过简单地去除落在“噪声空间”中的带噪信号向量分量来估计纯净信号。要将带噪信号向量空间分解为“信号子空间”和“噪声子空间”,可以采用线性代数中的正交矩阵分解技术,特别是奇异值分解(SVD)和特征值分解(EVD)。
二、源代码
clear;
% 调用MATLAB中含有噪声的数据文件 leleccum;
load leleccum;
index=1:3000;
x=leleccum(index);
N=8;
slength = length(x);
M=slength-100;
subplot(221);plot(x(1:M));
title('原始信号');
% 形成数据矩阵;
Signal=zeros(N,M);
for i=1:N
Signal(i,:)=x(i:M+i-1);
end
% 对数据矩阵作特征值分解;
[U, S, V]=svd(Signal);
d=diag(S(1:N,1:N));
subplot(222);stem(d);
title('特征值');
for i=1:N
if d(i)<mean(d)
d(i)=0;
end
end
stemp=S;
stemp(1:8,1:8)=diag(d);
Sf=U*stemp*V';
subplot(223);plot(Sf(1,:));
title('滤波之后的信号;阈值为特征值的平均值');
d=diag(S(1:N,1:N));
for i=1:N
if d(i)<=median(d)
d(i)=0;
end
end
% 声波信号降噪处理
function [TempData,Sample] =ReadData(FileName)
%数据读取,数据格式为特检院检测仪存储格式
% FileName='SR10-100021-SPL-5-4.dat';
fid=fopen(FileName,'rb');
%读取文件头
% 采样时间结构
Sample_Time = struct('year',{2009},'month',{1}, 'day',{22},'hour',{12},'minute',{12},'second',{12}) ;%定义采样时间格式
SensorAddr=fread(fid,1,'int32');%传感器地址
fseek(fid,4,'bof');%数据文件头长度为44字节
Sample_Time.year=fread(fid,1,'int32');
Sample_Time.month=fread(fid,1,'int32');
Sample_Time.day=fread(fid,1,'int32');
Sample_Time.hour=fread(fid,1,'int32');
Sample_Time.minute=fread(fid,1,'int32');
Sample_Time.second=fread(fid,1,'int32');
GroupOrder=fread(fid,1,'int8'); %读取的批次,0xff表示实时采样数据
startPack=fread(fid,1,'int8'); %起始包号,0~255
totalPack=fread(fid,1,'int8'); %总包数,0~255,0表示256包
location=fread(fid,1,'int8'); %/GPS定位信息
%设备配置信息结构 Device_Config
Device_Config.group=fread(fid,1,'int8'); %批号
Device_Config.year=fread(fid,1,'int8');%年
Device_Config.month=fread(fid,1,'int8');%月
Device_Config.day=fread(fid,1,'int8');%日
Device_Config.hour=fread(fid,1,'int8');%时
Device_Config.minute=fread(fid,1,'int8');%分
Device_Config.second=fread(fid,1,'int8');%秒
Device_Config.repeat=fread(fid,1,'int8');%次数
Device_Config.interval=fread(fid,1,'int8');%间隔
Device_Config.sampleLen=fread(fid,1,'int8');%采样长度
Device_Config.speed=fread(fid,1,'int8');%采样速度
Device_Config.gain=fread(fid,1,'int8');%增益
fseek(fid,44,'bof');%数据文件头长度为44字节
Count=261120;%数据长度
[TempData,Count]=fread(fid,'int16');
TempData=TempData*2.5/2048;
%得到数据长度
switch Device_Config.sampleLen
case {8}
Sample.Len='512KB';
case {9}
Sample.Len='512KB';
case {10}
Sample.Len='512B';
case {11}
Sample.Len='1024B';
end
% //采样速度
switch Device_Config.speed
case {1}
Sample.Frequency = '500k';
case {2}
Sample.Frequency = '250k';
case {3}
Sample.Frequency = '125k';
case {4}
Sample.Frequency = '62.5k';
end
% //增益
switch Device_Config.gain
case {0}
Sample.Gain = 1;
case {1}
Sample.Gain = 4;
case {2}
Sample.Gain = 16;
case {3}
Sample.Gain = 64;
case {4}
Sample.Gain = 10;
case {5}
Sample.Gain = 40;
case {6}
Sample.Gain = 160;
case {7}
Sample.Gain = 640;
case {8}
Sample.Gain = 100;
case {9}
Sample.Gain = 400;
case {10}
Sample.Gain = 1600;
case {11}
Sample.Gain = 6400;
end
三、运行结果
四、备注
版本:2014a
完整代码或代写加1564658423