📄 rice-rician distribution demo.htm
字号:
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
<!-- saved from url=(0079)http://www.mathworks.com/matlabcentral/files/14237/content/rician/ricedemo.html -->
<HTML xmlns:mwsh =
"http://www.mathworks.com/namespace/mcode/v1/syntaxhighlight.dtd"><HEAD><TITLE>Rice/Rician Distribution Demo</TITLE>
<META http-equiv=Content-Type content="text/html; charset=utf-8"><!--This HTML is auto-generated from an M-file.To make changes, update the M-file and republish this document. -->
<META content="MSHTML 6.00.2900.3059" name=GENERATOR>
<META content=2007-03-12 name=date>
<META content=ricedemo name=m-file>
<STYLE>BODY {
MARGIN: 10px; BACKGROUND-COLOR: white
}
H1 {
FONT-SIZE: x-large; COLOR: #990000
}
H2 {
FONT-SIZE: medium; COLOR: #990000
}
P {
WIDTH: 600px; ; WIDTH: expression(document.body.clientWidth > 620 ? "600px": "auto" ); max-width: 600px
}
H1 {
WIDTH: 600px; ; WIDTH: expression(document.body.clientWidth > 620 ? "600px": "auto" ); max-width: 600px
}
H2 {
WIDTH: 600px; ; WIDTH: expression(document.body.clientWidth > 620 ? "600px": "auto" ); max-width: 600px
}
DIV.content DIV {
WIDTH: 600px; ; WIDTH: expression(document.body.clientWidth > 620 ? "600px": "auto" ); max-width: 600px
}
PRE.codeinput {
PADDING-RIGHT: 10px; PADDING-LEFT: 10px; BACKGROUND: #eeeeee; PADDING-BOTTOM: 10px; PADDING-TOP: 10px
}
SPAN.keyword {
COLOR: #0000ff
}
SPAN.comment {
COLOR: #228b22
}
SPAN.string {
COLOR: #a020f0
}
SPAN.untermstring {
COLOR: #b20000
}
SPAN.syscmd {
COLOR: #b28c00
}
PRE.codeoutput {
PADDING-RIGHT: 10px; PADDING-LEFT: 10px; PADDING-BOTTOM: 10px; COLOR: #666666; PADDING-TOP: 10px
}
PRE.error {
COLOR: red
}
P.footer {
FONT-WEIGHT: lighter; FONT-SIZE: xx-small; COLOR: gray; FONT-STYLE: italic; TEXT-ALIGN: right
}
</STYLE>
</HEAD>
<BODY>
<DIV class=content>
<H1>Rice/Rician Distribution Demo</H1><INTRODUCTION>
<P>Demonstrate ricernd, ricepdf, and ricestat, in the context of simulating
Rician distributed noise for Magnetic Resonance Imaging magnitude data.
</P></INTRODUCTION>
<H2>Contents</H2>
<DIV>
<UL>
<LI><A
href="http://www.mathworks.com/matlabcentral/files/14237/content/rician/ricedemo.html#1">Example:
Rician noise vs (additive) Gaussian noise</A>
<LI><A
href="http://www.mathworks.com/matlabcentral/files/14237/content/rician/ricedemo.html#2">The
Rician PDF</A>
<LI><A
href="http://www.mathworks.com/matlabcentral/files/14237/content/rician/ricedemo.html#3">Approximations
to the Rician PDF</A>
<LI><A
href="http://www.mathworks.com/matlabcentral/files/14237/content/rician/ricedemo.html#6">A
pitfal: Adding approximate Rician noise</A> </LI></UL></DIV>
<H2>Example: Rician noise vs (additive) Gaussian noise<A name=1></A></H2>
<P>"Rician noise" depends on the data itself, it is not additive, so to "add"
Rician noise to data, what we really mean is make the data Rician distributed.
</P><PRE class=codeinput>clear, close(<SPAN class=string>'all'</SPAN>), clc
data = load(<SPAN class=string>'mri'</SPAN>); <SPAN class=comment>% MRI dataset (magnitude image) distributed with MATLAB</SPAN>
im = double(data.D(:,:,1,16)); <SPAN class=comment>% (particular slice through the brain)</SPAN>
s=5; <SPAN class=comment>% noise level (NB actual Rician stdev depends on signal, see ricestat)</SPAN>
im_g = im + 5 * randn(size(im)); <SPAN class=comment>% *Add* Gaussian noise</SPAN>
im_r = ricernd(im, s); <SPAN class=comment>% "Add" Rician noise (make Rician distributed)</SPAN>
<SPAN class=comment>% Note that the Gaussian noise has made some data go negative.</SPAN>
min_o = round(min(im(:))); max_o = round(max(im(:)));
min_g = round(min(im_g(:))); max_g = round(max(im_g(:)));
min_r = round(min(im_r(:))); max_r = round(max(im_r(:)));
<SPAN class=comment>% Show each image with the same color scaling limits</SPAN>
clim = [min_g max(max_g, max_r)];
figure(<SPAN class=string>'Position'</SPAN>, [30 500 800 300]); colormap(data.map);
subplot(1,3,1); imagesc(im, clim); axis <SPAN class=string>image</SPAN>;
title(<SPAN class=string>'Original'</SPAN>);
xlabel([<SPAN class=string>'Range: ('</SPAN> num2str(min_o) <SPAN class=string>', '</SPAN> num2str(max_o) <SPAN class=string>')'</SPAN>])
subplot(1,3,2); imagesc(im_g, clim); axis <SPAN class=string>image</SPAN>;
title(<SPAN class=string>'Gaussian noise'</SPAN>);
xlabel([<SPAN class=string>'Range: ('</SPAN> num2str(min_g) <SPAN class=string>', '</SPAN> num2str(max_g) <SPAN class=string>')'</SPAN>])
subplot(1,3,3); imagesc(im_r, clim); axis <SPAN class=string>image</SPAN>;
title(<SPAN class=string>'Rician noise'</SPAN>)
xlabel([<SPAN class=string>'Range: ('</SPAN> num2str(min_r) <SPAN class=string>', '</SPAN> num2str(max_r) <SPAN class=string>')'</SPAN>])
</PRE><IMG hspace=5 src="Rice-Rician Distribution Demo.files/ricedemo_01.png"
vspace=5>
<H2>The Rician PDF<A name=2></A></H2><PRE class=codeinput>x = linspace(0, 8, 100);
figure;
subplot(3, 1, 1)
plot(x, ricepdf(x, 0, 1), x, ricepdf(x, 1, 1),<SPAN class=keyword>...</SPAN>
x, ricepdf(x, 2, 1), x, ricepdf(x, 4, 1))
title(<SPAN class=string>'Rice PDF with s = 1'</SPAN>)
legend(<SPAN class=string>'v=0'</SPAN>, <SPAN class=string>'v=1'</SPAN>, <SPAN class=string>'v=2'</SPAN>, <SPAN class=string>'v=4'</SPAN>)
subplot(3,1,2)
plot(x, ricepdf(x, 1, 0.25), x, ricepdf(x, 1, 0.50),<SPAN class=keyword>...</SPAN>
x, ricepdf(x, 1, 1.00), x, ricepdf(x, 1, 2.00))
title(<SPAN class=string>'Rice PDF with v = 1'</SPAN>)
legend(<SPAN class=string>'s=0.25'</SPAN>, <SPAN class=string>'s=0.50'</SPAN>, <SPAN class=string>'s=1.00'</SPAN>, <SPAN class=string>'s=2.00'</SPAN>)
subplot(3,1,3)
plot(x, ricepdf(x, 0.5, 0.25), x, ricepdf(x, 1, 0.50),<SPAN class=keyword>...</SPAN>
x, ricepdf(x, 2.0, 1.00), x, ricepdf(x, 4, 2.00))
title(<SPAN class=string>'Rice PDF with v/s = 2'</SPAN>)
legend(<SPAN class=string>'s=0.25'</SPAN>, <SPAN class=string>'s=0.50'</SPAN>, <SPAN class=string>'s=1.00'</SPAN>, <SPAN class=string>'s=2.00'</SPAN>)
</PRE><IMG hspace=5 src="Rice-Rician Distribution Demo.files/ricedemo_02.png"
vspace=5>
<H2>Approximations to the Rician PDF<A name=3></A></H2>
<P>At high SNR (e.g. v/s > 3), Rician data is approximately Gaussian, though
note that v is not the mean of the Rician distribution, nor of the Gaussian
approximation (though for very high SNR it tends towards it). </P><PRE class=codeinput>v = 1.75; s = 0.5; x = linspace(0, 4, 100);
mu1 = v;
gpdf1 = (2*pi*s^2)^(-0.5) * exp(-0.5*(x-mu1).^2/s^2);
mu2 = ricestat(v, s); <SPAN class=comment>% (mu2 is approximately sqrt(v^2 + s^2))</SPAN>
gpdf2 = (2*pi*s^2)^(-0.5) * exp(-0.5*(x-mu2).^2/s^2);
figure; plot(x, ricepdf(x, v, s), x, gpdf1, x, gpdf2);
title(<SPAN class=string>'High SNR'</SPAN>);
legend(<SPAN class=string>'Rician'</SPAN>, <SPAN class=string>'Gaussian, naive mean'</SPAN>, <SPAN class=string>'Gaussian, corrected mean'</SPAN>)
</PRE><IMG hspace=5 src="Rice-Rician Distribution Demo.files/ricedemo_03.png"
vspace=5>
<P>At very low SNR, Rician data is approximately Rayleigh distributed (with no
signal, the distributions are exactly equivalent)</P><PRE class=codeinput>v = 0.25; s = 0.5; x = linspace(0, 3, 100);
raypdf = (x / s^2) .* exp(-x.^2 / (2*s^2));
figure; plot(x, ricepdf(x, v, s), x, raypdf);
title(<SPAN class=string>'Low SNR'</SPAN>); legend(<SPAN class=string>'Rician'</SPAN>, <SPAN class=string>'Rayleigh'</SPAN>)
</PRE><IMG hspace=5 src="Rice-Rician Distribution Demo.files/ricedemo_04.png"
vspace=5>
<P>At low-medium SNR, neither Gaussian nor Rayleigh is a great approximation</P><PRE class=codeinput>v = 0.75; s = 0.5; x = linspace(-1, 3, 100);
mu = ricestat(v, s);
gpdf = (2*pi*s^2)^(-0.5) * exp(-0.5*(x-mu).^2/s^2);
raypdf = (x / s^2) .* exp(-x.^2 / (2*s^2)); raypdf(x <= 0) = 0;
figure; plot(x, ricepdf(x, v, s), x, gpdf, x, raypdf);
title(<SPAN class=string>'Medium SNR'</SPAN>); legend(<SPAN class=string>'Rician'</SPAN>, <SPAN class=string>'Gaussian'</SPAN>, <SPAN class=string>'Rayleigh'</SPAN>)
</PRE><IMG hspace=5 src="Rice-Rician Distribution Demo.files/ricedemo_05.png"
vspace=5>
<H2>A pitfal: Adding approximate Rician noise<A name=6></A></H2>
<P>Since the Rician distribution with zero signal is equivalent to the Rayleigh,
and with high SNR is approximated by a Gaussian, it is tempting to add Rayleigh
or Gaussian noise (depending on SNR) to existing data, to simulate Rician
distributed data (aka "adding Rician noise"). </P><PRE class=codeinput>v = 75; s = 50; <SPAN class=comment>% medium SNR => Add Rayleigh noise</SPAN>
<SPAN class=comment>% Generate Rayleigh samples</SPAN>
N = 10000;
gsamp = s * randn(2, N);
rsamp = sqrt(sum(gsamp .^ 2));
</PRE>
<P>Compare data with added Rayleigh samples to expected Rician distribution</P><PRE class=codeinput>r = v + rsamp; <SPAN class=comment>% Add Rayleigh noise to (constant) signal</SPAN>
c = linspace(0, 250, 20);
w = c(2); <SPAN class=comment>% histogram bin-width</SPAN>
h = histc(r, c);
figure; hold <SPAN class=string>on</SPAN>; bar(c, h, <SPAN class=string>'histc'</SPAN>);
xl = xlim; x = linspace(xl(1), xl(2), 100);
yl = ylim; <SPAN class=comment>% (for second plot)</SPAN>
plot(x, N*w*ricepdf(x, v, s), <SPAN class=string>'r'</SPAN>);
title(<SPAN class=string>'Histogram of data, vs desired Rician PDF'</SPAN>)
legend(<SPAN class=string>'Data with added Rayleigh noise'</SPAN>, <SPAN class=string>'Rician distribution'</SPAN>)
</PRE><IMG hspace=5 src="Rice-Rician Distribution Demo.files/ricedemo_06.png"
vspace=5>
<P>Compare Rician samples to expected Rician distribution</P><PRE class=codeinput>r = ricernd(v*ones(1, N), s);
c = linspace(0, 250, 20);
w = c(2); <SPAN class=comment>% histogram bin-width</SPAN>
h = histc(r, c);
figure; hold <SPAN class=string>on</SPAN>; bar(c, h, <SPAN class=string>'histc'</SPAN>);
xl = xlim; x = linspace(xl(1), xl(2), 100);
ylim(yl);
plot(x, N*w*ricepdf(x, v, s), <SPAN class=string>'r'</SPAN>);
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -