matlab_频率选择性rician衰落信道

%% Rice/Rician Distribution Demo % Demonstrate ricernd, ricepdf, and ricestat, in the context of simulating % Rician distributed noise for Magnetic Resonance Imaging magnitude data. %% Example: Rician noise vs (additive) Gaussian noise % "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. clear, close('all'), clc data = load('mri'); % MRI dataset (magnitude image) distributed with MATLAB im = double(data.D(:,:,1,16)); % (particular slice through the brain) s=5; % noise level (NB actual Rician stdev depends on signal, see ricestat) im_g = im + 5 * randn(size(im)); % *Add* Gaussian noise im_r = ricernd(im, s); % "Add" Rician noise (make Rician distributed) % Note that the Gaussian noise has made some data go negative. 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(:))); % Show each image with the same color scaling limits clim = [min_g max(max_g, max_r)]; figure('Position', [30 500 800 300]); colormap(data.map); subplot(1,3,1); imagesc(im, clim); axis image; title('Original'); xlabel(['Range: (' num2str(min_o) ', ' num2str(max_o) ')']) subplot(1,3,2); imagesc(im_g, clim); axis image; title('Gaussian noise'); xlabel(['Range: (' num2str(min_g) ', ' num2str(max_g) ')']) subplot(1,3,3); imagesc(im_r, clim); axis image; title('Rician noise') xlabel(['Range: (' num2str(min_r) ', ' num2str(max_r) ')']) %% The Rician PDF x = linspace(0, 8, 100); figure; subplot(3, 1, 1) plot(x, ricepdf(x, 0, 1), x, ricepdf(x, 1, 1),... x, ricepdf(x, 2, 1), x, ricepdf(x, 4, 1)) title('Rice PDF with s = 1') legend('v=0', 'v=1', 'v=2', 'v=4') subplot(3,1,2) plot(x, ricepdf(x, 1, 0.25), x, ricepdf(x, 1, 0.50),... x, ricepdf(x, 1, 1.00), x, ricepdf(x, 1, 2.00)) title('Rice PDF with v = 1') legend('s=0.25', 's=0.50', 's=1.00', 's=2.00') subplot(3,1,3) plot(x, ricepdf(x, 0.5, 0.25), x, ricepdf(x, 1, 0.50),... x, ricepdf(x, 2.0, 1.00), x, ricepdf(x, 4, 2.00)) title('Rice PDF with v/s = 2') legend('s=0.25', 's=0.50', 's=1.00', 's=2.00') %% Approximations to the Rician PDF % 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). 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); % (mu2 is approximately sqrt(v^2 + s^2)) 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('High SNR'); legend('Rician', 'Gaussian, naive mean', 'Gaussian, corrected mean') %% % At very low SNR, Rician data is approximately Rayleigh distributed % (with no signal, the distributions are exactly equivalent) 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('Low SNR'); legend('Rician', 'Rayleigh') %% % At low-medium SNR, neither Gaussian nor Rayleigh is a great approximation 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('Medium SNR'); legend('Rician', 'Gaussian', 'Rayleigh') %% A pitfal: Adding approximate Rician noise % 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"). v = 75; s = 50; % medium SNR => Add Rayleigh noise % Generate Rayleigh samples N = 10000; gsamp = s * randn(2, N); rsamp = sqrt(sum(gsamp .^ 2)); %% % Compare data with added Rayleigh samples to expected Rician distribution r = v + rsamp; % Add Rayleigh noise to (constant) signal c = linspace(0, 250, 20); w = c(2); % histogram bin-width h = histc(r, c); figure; hold on; bar(c, h, 'histc'); xl = xlim; x = linspace(xl(1), xl(2), 100); yl = ylim; % (for second plot) plot(x, N*w*ricepdf(x, v, s), 'r'); title('Histogram of data, vs desired Rician PDF') legend('Data with added Rayleigh noise', 'Rician distribution') %% % Compare Rician samples to expected Rician distribution r = ricernd(v*ones(1, N), s); c = linspace(0, 250, 20); w = c(2); % histogram bin-width h = histc(r, c); figure; hold on; bar(c, h, 'histc'); xl = xlim; x = linspace(xl(1), xl(2), 100); ylim(yl); plot(x, N*w*ricepdf(x, v, s), 'r'); title('Histogram of data, vs desired Rician PDF') legend('Rician sampled data', 'Rician distribution') %% % We see that the additive Rayleigh result is very poor (with medium SNR). % % With high noise levels the appearance of results can differ dramatically: s = 0.15 * max_o; % Approximate additive Rayleigh/Gaussian noise: im_r = im; for i=1:numel(im_r) if im_r(i) <= 3*s im_r(i) = im_r(i) + s*norm(randn(2,1)); % Add Rayleigh noise else im_r(i) = im_r(i) + s*randn; % Add Gaussian noise end end % "Add" Rician noise (make Rician distributed): im_R = ricernd(im, s); % Show each image with the same color scaling limits min_r = round(min(im_r(:))); max_r = round(max(im_r(:))); min_R = round(min(im_R(:))); max_R = round(max(im_R(:))); clim = [min_r max(max_r, max_R)]; figure('Position', [30 500 800 300]); colormap(data.map); subplot(1,3,1); imagesc(im, clim); axis image; title('Original'); xlabel(['Range: (' num2str(min_o) ', ' num2str(max_o) ')']) subplot(1,3,2); imagesc(im_r, clim); axis image; title('Rayleigh/Gaussian noise') xlabel(['Range: (' num2str(min_r) ', ' num2str(max_r) ')']) subplot(1,3,3); imagesc(im_R, clim); axis image; title('Rician data'); xlabel(['Range: (' num2str(min_R) ', ' num2str(max_R) ')']) %% % Copyright 2007 Ged Ridgway (Ged at cantab dot net)