matlab对图像去噪后怎么求信噪比

Posted 2023-05-17

一般用峰值信噪比PSNR来衡量滤波降噪效果
其定义为:

对于8位图像,MAX=255,
其中MSE为均方误差:

I为原图像,K为滤波后图像 参考技术A matlab对图像去噪后求信噪比
处理的程序
X=imread('1.tif')
X=im2double(X)
nbc=size(map,1)
x=imnoise(X,'gaussian',0,0.005)
wname='bior3.3'
[c,s]=wavedec2(x,3,wname)
n=[1,2,3];p=[150,120,60]
nc1=wthcoef2('d',c,s,n,p,'s')
nc2=wthcoef2('d',c,s,n,p,'h')
rx1=waverec2(nc1,s,wname)
rx2=waverec2(nc2,s,wname)
[a,b]=psnr_mse(rx1,X)
colormap(gray(nbc))
subplot(221)
image(wcodemat(X,nbc))
title('原图')
subplot(222)
image(wcodemat(x,nbc))
title('含方差为0.005的高斯噪声')
subplot(223)
image(wcodemat(rx1,nbc))
title('小波基bior3.3软阈值去噪')
subplot(224)
image(wcodemat(rx2,nbc))
title('小波基bior3.3硬阈值去噪')
求psnr的程序
function
[mse,psnr]
=
psnr_mse(image,image_prime)
%
convert
to
doubles
image=double(image)
[M,N]=size(image)
image_prime=double(image_prime)
%
avoid
divide
by
zero
nastiness
if
((sum(sum((image-image_prime).^2)))
==
0)
error('Input
vectors
must
not
be
identical')
else
mse=sum(sum((image-image_prime).^2))/(M*N);
psnr=10*log10(255^2/mse)
end
return
编程语言（programming
language），又称程序设计语言（program
design
language，简称PDL），是一组用来定义计算机程序的语法规则。它是一种被标准化的交流技巧，用来向计算机发出指令。一种计算机语言让程序员能够准确地定义计算机所需要使用的数据，并精确地定义在不同情况下所应当采取的行动。计算机语言的种类非常的多，总的来说可以分成机器语言，汇编语言，高级语言三大类。
目前通用的编程语言有两种形式：汇编语言和高级语言。

语音去噪基于matlab GUI谱减法+维纳滤波语音去噪（带面板+信噪比）含Matlab源码 1661期

一、简介

1 维纳滤波法
维纳滤波法(wiener filter)也是一个比较经典的传统做法，它的本质是估计出一个线性滤波器，也就是一个向量，这个滤波器会对不同的频段进行不同程度的抑制，其保真效果会比谱减法要好一些。

我们这里不会讲详细的推导过程，只讲其大致思想。因为这么大功夫推导出来，还是有很多不能解决的问题，还不如深度学习train一发。想看详细推导了可以去看知乎的卡尔曼滤波器详解——从零开始(3) Kalman Filter from Zero这篇，于泓-语音增强-维纳滤波这个视频讲的更偏向于应用，都很棒。

还有就是这里讲的是smoothing的问题，即根据未来的信号，过去的信号以及现在的信号来推测出现在的干净信号。除此之外，还有prediction和filtering的问题，prediction是指根据过去的和现在的信号，预测未来的干净信号；filtering的问题是指根据过去和现在的信号，推测现在的干净信号。所以这里讲的方法没法应用于实时语音去噪，只能在拿到整段信号之后，对这段信号进行去噪。

维纳滤波器的设计准则为使得干净信号x(n)x(n)和估计的干净信号x ^ ( n ) 之间的差值越小越好，即计算一个最小均方差

2 谱减法
应该是最早被用于语音去噪的方法，它的思想非常简单，就是通过估计出噪声，并在频域里将其幅值剪掉，再还原，就结束了。为了表示方便，我们先假设纯净的声音为x(n)，原始声音为y(n)，噪声为e(x)，就有y(n)=x(n)+e(n)
这里只有y(n)是我们有的，其他x(n)和e(n)都还不知道，目的是把x(n)给求出来。
noisereduce中的stationary的方法就是用谱减法去做的，效果还是不错的，不过也只能应对于stationary noise。
我们按谱减法步骤来说明一下整个过程。
（1）截取头部一小段语音作为噪声
e(m)=y(n)[:m]
其中e表示噪声信号，y表示原始信号，m和n表示sample的数量。
我们认为stationary noise是一直存在于背景当中的声音，而人声一般在开头的几十毫秒是没有的，所以就默认取前面一小段作为噪声。不过当无法确定噪声的时候，把整段声音都作为噪声也是可以的，noisereduce就是这么做的。
（2）分别计算原始音频和噪声的STFT，Y(ω)和E(ω)。
（3）根据噪声的频谱幅值，对原始音频的频谱幅值进行谱减。

这样做不好的地方就是会有很多坑坑洼洼的噪声频率残留，这个现象也被称为是音乐噪声。实际操作过会发现这种方法减了和没减差不多。因此有人提出了过减法，就是宁可错杀一千不能放过一个的做法。

其中，α∈[0,+∞]是过减因子，β∈[0,1]是谱下限参数，用来控制残留多少的噪声。这样减出来噪声会明显少了很多，但声音也会随着α \\alphaα的增大而逐渐失真。
noisereduce中的具体实现略有不同，它过减用∣E(ω)∣的方差来控制，一般是1.5倍或者1.0倍的方差。代码片段如下所示

self.mean_freq_noise = np.mean(noise_stft_db, axis=1)
self.std_freq_noise = np.std(noise_stft_db, axis=1)
self.noise_thresh = self.mean_freq_noise + self.std_freq_noise * self.n_std_thresh_stationary

小于noise_thresh的幅值会置0，其余的保留。n_std_thresh_stationary为0时，就是没有过减的式(2−2)。
（4）对∣ X ^ ( ω ) ∣做平滑处理，使得声音失真没那么严重。
noisereduce中使用的scipy.signal.fftconvolve来实现这一过程。
（5）结合原始音频的相位，还原谱减后的音频。这就是个反向STFT的过程。

二、部分源代码

% algorithm_flag = 0;
function varargout = gui(varargin)
% GUI MATLAB code for gui.fig
%      GUI, by itself, creates a new GUI or raises the existing
%      singleton*.
%
%      H = GUI returns the handle to a new GUI or the handle to
%      the existing singleton*.
%
%      GUI('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in GUI.M with the given input arguments.
%
%      GUI('Property','Value',...) creates a new GUI or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before gui_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to gui_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help gui

% Last Modified by GUIDE v2.5 06-Apr-2021 19:48:53

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @gui_OpeningFcn, ...
                   'gui_OutputFcn',  @gui_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin1)
    gui_State.gui_Callback = str2func(varargin1);
end

if nargout
    [varargout1:nargout] = gui_mainfcn(gui_State, varargin:);
else
    gui_mainfcn(gui_State, varargin:);
end
% End initialization code - DO NOT EDIT


% --- Executes just before gui is made visible.
function gui_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to gui (see VARARGIN)

% Choose default command line output for gui
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes gui wait for user response (see UIRESUME)
% uiwait(handles.figure1);


% --- Outputs from this function are returned to the command line.
function varargout = gui_OutputFcn(hObject, eventdata, handles) 
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout1 = handles.output;


%算法1
function pushbutton1_Callback(hObject, eventdata, handles)
set(handles.text6, 'String','谱减法');
%谱减法语音降噪
global filename;
[xx, fs] = audioread(filename);           % 读入数据文件  样本数据 xx  该数据的采样率 Fs。
xx=xx-mean(xx);                         % 消除直流分量
global x;
x=xx/max(abs(xx));                      % 幅值归一化   abs 数据绝对值

IS=0.25;                                % 设置前导无话段长度
wlen=200;                               % 设置帧长为25ms
inc=80;                                 % 设置帧移为10ms
global SNR;
%SNR=5;                                  % 设置信噪比SNR
N=length(x);                            % 信号长度

global signal
signal=awgn(x,SNR,'measured','db');               % 叠加噪声

global snr1;
global snr2;
snr1=SNR_Calc(x,signal);            % 计算初始信噪比
NIS=fix((IS*fs-wlen)/inc +1);           % 求前导无话段帧数

a=4; b=0.001;                           % 设置参数a和b
global output;
output=SpectralSub(signal,wlen,inc,NIS,a,b);% 谱减
snr2=SNR_Calc(x,output);            % 计算谱减后的信噪比
snr=snr2-snr1;
fprintf('snr1=%5.4f   snr2=%5.4f   snr=%5.4f\\n',snr1,snr2,snr);       
% %显示图片
% handles.x = x;
% axes(handles.axes1);
% plot(handles.x);    
% handles.signal = signal;
% axes(handles.axes2);
% plot(handles.signal);
% handles.output = output;
% axes(handles.axes3);
% plot(handles.output);


%算法2
function pushbutton2_Callback(hObject, eventdata, handles)
set(handles.text6, 'String','维纳滤波');
%基于先验信噪比的维纳滤波算法语音降噪
global filename;
[xx, fs] = audioread(filename);           % 读入数据文件
xx=xx-mean(xx);                         % 消除直流分量
global x;
x=xx/max(abs(xx));                      % 幅值归一化

IS=0.25;                                % 设置前导无话段长度
wlen=200;                               % 设置帧长为25ms
inc=80;                                 % 设置帧移为10ms
global SNR;
%SNR=5;                                  % 设置信噪比SNR
NIS=fix((IS*fs-wlen)/inc +1);           % 求前导无话段帧数
alpha=0.95;
global signal;
signal=awgn(x,SNR,'measured','db');               % 叠加噪声

global output;
output=Weina_Im(x,wlen,inc,NIS,alpha) ;
output=output/max(abs(output));
len=min(length(output),length(x));
x=x(1:len);
signal=signal(1:len);
output=output(1:len);
global snr1;
global snr2;
snr1=SNR_Calc(x,signal);            % 计算初始信噪比
snr2=SNR_Calc(x,output);            % 计算降噪后的信噪比
snr=snr2-snr1;
fprintf('snr1=%5.4f   snr2=%5.4f   snr=%5.4f\\n',snr1,snr2,snr);
% %显示图片
% handles.x = x;
% axes(handles.axes1);
% plot(handles.x);    
% handles.signal = signal;
% axes(handles.axes2);
% plot(handles.signal);
% handles.output = output;
% axes(handles.axes3);
% plot(handles.output);

三、运行结果

四、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1]韩纪庆,张磊,郑铁然.语音信号处理（第3版）[M].清华大学出版社，2019.
[2]柳若边.深度学习:语音识别技术实践[M].清华大学出版社，2019.