【数字信号去噪】基于matlab奇异值分解（SVD）数字信号降噪【含Matlab源码 1020期】

原创

佛怒火莲 2022-04-07 16:20:52 博主文章分类：matlab完整代码 ©著作权

©著作权归作者所有：来自51CTO博客作者佛怒火莲的原创作品，请联系作者获取转载授权，否则将追究法律责任

二、简介

基于SVD（奇异值分解）的去噪声技术属于子空间算法的一种。简单的来说我们希望将带噪信号向量空间分解为分别由纯净信号主导和噪声信号主导的两个子空间，然后通过简单地去除落在“噪声空间”中的带噪信号向量分量来估计纯净信号。要将带噪信号向量空间分解为“信号子空间”和“噪声子空间”，可以采用线性代数中的正交矩阵分解技术，特别是奇异值分解（SVD）和特征值分解（EVD）。

【数字信号去噪】基于matlab奇异值分解（SVD）数字信号降噪【含Matlab源码 1020期】_数据

【数字信号去噪】基于matlab奇异值分解（SVD）数字信号降噪【含Matlab源码 1020期】_matlab_02

【数字信号去噪】基于matlab奇异值分解（SVD）数字信号降噪【含Matlab源码 1020期】_matlab_03

【数字信号去噪】基于matlab奇异值分解（SVD）数字信号降噪【含Matlab源码 1020期】_matlab_04

【数字信号去噪】基于matlab奇异值分解（SVD）数字信号降噪【含Matlab源码 1020期】_自动驾驶_05

【数字信号去噪】基于matlab奇异值分解（SVD）数字信号降噪【含Matlab源码 1020期】_特征值_06

三、部分源代码

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