% file RMS.m
%
% Program to find the RMS value for the random array to be used later
%
% by X. R. Wang, 12th May 2000
global SIZE N sigma2 % global variables
%SIZE = 128; % uncomment to test this program
SIZE2 = SIZE^2;
% load the data array. four such arrays are available
% uncomment the appropriate line to load the file.
% load straight % straight vertical line
% load diagonal % diagonal line
% load several % several lines with different gradients
% load square % square
% straight % uncomment to generate the array here
% change between straight, diagonal, several and square
% to generate different noise
% turn the array into an image to be placed in figure window 1
figure(1)
colormap(jet) % using jet colour scheme for the image
image(N*50) % enhance every element of the array 50 times stronger
zoom on
sumN = sum(sum(N)); % find the sum of all the elements in the matrix N
sqN = N.^2; % square every element of matrix N
sum_sq = sum(sum(sqN)); % find the sum of the square of N(i,j)
avg = sumN / SIZE2; % find the mean of the array
varience = sum_sq / SIZE2 - avg*avg; % find the value of varience
Wiener filter
The final program is one with the actual wiener filtering and it is as follows:
% wiener.m
%
% Program to wiener filtering signal from an array of noise
% need RMS.m and either data array in a file or the code to generate the array
%
% Written for PHYS3301 Scientific Computing
% by X. Rosalind Wang, 19th May, 2000
global SIZE N sigma % global variables
SIZE = 128; % size of the noise array
BOX = 5; % size of the box, need to be an odd number
BOX2 = (BOX-1)/2;
START = BOX2 + 1;
END = SIZE - START;
BBOX = BOX*BOX;
RMS % call function RMS.m to load array and find the RMS value
n=1; m=1; % counter
B=zeros(SIZE); % array of zeros (no signal or noise),
% same size as the original array
for i = START:END
for j = START:END
% Copy the content in each box into a different array, and find its properties
M = N(i-BOX2:i+BOX2,j-BOX2:j+BOX2);
a = sum(sum(M));
M2 = M.^2;
b = sum(sum(M2));
mu = a/BBOX; % find local mean
sig = b/BBOX - mu*mu; % find local variance
sig = max(sig,10e-8);
% filtering
A(n:n+BOX-1,m:m+BOX-1) = mu+((sig-varience)/sig)*(N(i-BOX2:i+BOX2,j-BOX2:j+BOX2)-mu);
B(n:n+BOX-1,m:m+BOX-1) = (B(n:n+BOX-1,m:m+BOX-1) + A(n:n+BOX-1,m:m+BOX-1)) ./2;
m = m+1 % move box to the right by one
end
n = n+1; % go to next row
m = 1; % start from column one
end
% print out the filtered signal on figure(2)
figure(2)
colormap(jet)
image(B*50)
zoom on
%I have chose the jet colour to display the data array because it has very clear
%and bright colour contrast from red to orange to yellow to light blue and finally to dark blue.
%This colour contrast shows very clearly the signal from the noise in the filtered data array.