【數字信号去噪】奇異值分解（SVD）數字信号降噪【含Matlab源碼 1020期】

一、簡介

基于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