## Pages

### Matlab Sample Codes

1. Performance Analysis of BPSK in the presence of AWGN Noise and Rayleigh Fading Channel

x_bits= randi([0 1],1,1000000);   %   Input  bits
hq = modem.pskmod('M', 2, 'PhaseOffset', 0, 'SymbolOrder', 'Gray', 'InputType', 'Bit');
%Modulation Object for      BPSK
x=modulate(hq, x_bits);                  %    BPSK symbols
LTx=length(x);                            %   length of symbol vector
Eb_No=[0:1:10];                          %   energy per bit talen in (dB) which is a conventional way
LE=length(Eb_No);

% Eb/No;                                 %  This line does not make ay sense  makes no sense

%Es_No=Eb_No+10*log10(1/2); appropriate relation is Es/N=(Eb/N)*m where m
%is number of bits per symbol tatat depends on
%modultion order. For e.g. for BPSK these two
%are same , for QPSK or QAM Es/N=(Eb/N)*2. When
%we convert in (dB) we may write as:

%   Es_No=Eb_No+10*log10(m);

Es_No=Eb_No;   %  BPSK

% Es_No (dB) = 10 log10(Symbol_POwer/Noise_Power)
%Symbol_Power/Noise_Power=10^(Es_No./10)
%Usually symbol power is normalized i.e. Symbol_Power=1 in M-PSK or others
No=1./(10.^(Es_No./10)); % Noise power or Noise variance (sigma square)

ber=[];
for i=1:LE

% m+s*randn(1,N) genrates N random numbers with mean as (m) and standard
% deviation as (s). see matlab help for randn
% As we lnow white noise have zero mean(m=0) and variance as No (s=sqrt(No);

% is you use  Es_No=Eb_No+10*log10(m); then appropriate expression shpoul
% be,   nx=sqrt(No(i)*m)*randn(1,LTx);

nx=sqrt(No(i))*randn(1,LTx);  % For BPSK
% y=x+nx(1:length(x)); it is better to generate noise vector equal to
% symbol vector. Here in this case these are same so we can write:

y=x+nx;

hq = modem.pskdemod('M', 2, 'PhaseOffset', 0,'SymbolOrder','Gray', 'OutputType', 'Bit','DecisionType', 'hard decision');
y_bits = demodulate(hq, y);

x_bits=x_bits(:);
y_bits=y_bits(:);

s=0;
for j=1:length(x_bits)
if x_bits(j) ~= y_bits(j)
s=s+1;
end
end

br=s/length(x_bits);
ber=[ber br];
end
semilogy(Es_No,ber);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Rayleigh

ber=[];
for i=1:LE

nx=sqrt(No(i))*randn(1,LTx);  % For BPSK

% rayleigh_chn=raylrnd(1,1,LTx);   %   works in MATLAB 2013 , dontt know
rayleigh_chn=abs(randn(1,LTx)+1j*randn(1,LTx));

y=x.*rayleigh_chn+nx;

hq = modem.pskdemod('M', 2, 'PhaseOffset', 0,'SymbolOrder','Gray', 'OutputType', 'Bit','DecisionType', 'hard decision');
y_bits = demodulate(hq, y);

x_bits=x_bits(:);
y_bits=y_bits(:);

s=0;
for j=1:length(x_bits)
if x_bits(j) ~= y_bits(j)
s=s+1;
end
end

br=s/length(x_bits);
ber=[ber br];
end
hold on;
semilogy(Es_No,ber,'r');

2. Matlab Code for Wiener Filtering of an Image in Frequency Domain

function [fgzwnr] = wienerindft(fgz)

% fgz  input image;

v=std2(fgz);   %standard deviation of input image
v0=v^2

[M,N]=size(fgz);
F = fft2(fgz);
%Fmag = abs(F);        %  magnitude
%Fmag = abs(F);        % magnitude
Fmag = abs(F/sqrt(M*N));        %  normalized magnitude

Fmag1 = modified1(Fmag, v0);

fzero = find(Fmag==0); Fmag(fzero)=1; Fmag1(fzero)=0;
F = F.*Fmag1./Fmag;

% inverse FFT transform
fgzwnr = real(ifft2(F));

end

function [out]=modified1(in,v0)
S=size(in);
%in=double(in);
%v0=25;
W= [3,5,7,9];
%//////////////////////////////////////////////////////////////////////////
N1=ones(W(1),W(1));
N2=ones(W(2),W(2));
N3=ones(W(3),W(3));
N4=ones(W(4),W(4));
%//////////////////////////////////////////////////////////////////////////

sigma1 = conv2(in.*in,N1,'same')/(W(1)*W(1))- v0;
N2=ones(W(2),W(2));
sigma2 = conv2(in.*in,N2,'same')/(W(2)*W(2))- v0;
N3=ones(W(3),W(3));
sigma3 = conv2(in.*in,N3,'same')/(W(3)*W(3))- v0;
N4=ones(W(4),W(4));
sigma4 = conv2(in.*in,N4,'same')/(W(4)*W(4))- v0;

%//////////////////////////////////////////////////////////////////////////

for i= 1:S(1)
for j=1:S(2)
sigma1(i,j)= max(0,sigma1(i,j));
sigma2(i,j)= max(0,sigma2(i,j));
sigma3(i,j)= max(0,sigma3(i,j));
sigma4(i,j)= max(0,sigma4(i,j));
a= [sigma1(i,j),sigma2(i,j),sigma3(i,j),sigma4(i,j)];
sigma(i,j)= min(a);
out(i,j)= in(i,j)* sigma(i,j)/(sigma(i,j)+v0);
end
end
end

3. Singular Value Decomposition (SVD) Based Image Compression

function []= SVD_C()

k=10;
% determine size of host image
Mc=size(hostimage,1);            %Height
Nc=size(hostimage,2);            %Width

[Uh,Sh,Vh]= svd(double(hostimage));

VhT=transpose(Vh(1:Nc,1:k));
compressed_image= Uh(1:Mc,1:k)*Sh(1:k,1:k)*VhT;

%PSNR

[PSNR_SVD,MSE_SVD]= psnr(hostimage,compressed_image)

imwrite(compressed_image,'lena_svd.jpg','jpg');

figure(1)
imshow(hostimage,[]);
title('Host Image');
figure(2)
imshow(compressed_image,[])
title('Compressed Image');

end

% Function to calculate psnr used above
function [p mse]= psnr(original,reconstructed)

Mc= size(original,1);
Nc= size(original,2);
a1= original(1:Mc,1:Nc);
a2= reconstructed(1:Mc,1:Nc);

error=double(a1)-double(a2);

sum=0;
for (i=1:Mc)
for(j=1:Nc)
sum = sum + ((error(i,j))^2);
end
end

mse = sum/(Mc*Nc);
p = 10*log(double((255*255)/mse));

end

4. Image Compression Using Discrete Cosine Transform (DCT)

function []= DCT_C()
blocksize=8; % set the size of the block
quant_multiple = 1;

% determine size of host image
Mc=size(hostimage,1); %Height
Nc=size(hostimage,2); %Width

blockcount = (Mc*Nc)/8^2;

DCT_quantizer = [ 16 11 10 16 24 40 51 61;
12 12 14 19 26 58 60 55;
14 13 16 24 40 57 69 56;
14 17 22 29 51 87 80 62;
18 22 37 56 68 109 103 77;
24 35 55 64 81 104 113 92;
49 64 78 87 103 121 120 101;
72 92 95 98 112 100 103 99 ];

% process the image in blocks

x=1;
y=1;
y1=1;
x1=1;

for (i = 1:blockcount)

% transform block using DCT
dct_block=dct2(hostimage(y:y+blocksize-1,x:x+blocksize-1));

% Quantization
dct_block = floor (dct_block ./ (DCT_quantizer(1:blocksize, 1:blocksize) * quant_multiple) + 0.5);

% transform block back into spatial domain
compressed_image(y:y+blocksize-1,x:x+blocksize-1)=idct2(dct_block);

% move on to next block. At end of row move to next row
if (x+blocksize) > Nc
x=1;
y=y+blocksize;

else
x=x+blocksize;
end

end

%PSNR
[PSNR_DCT,MSE_DCT]= psnr(hostimage,compressed_image)  % psnr function code written below

imwrite(compressed_image,'lena_dct.jpg','jpg');

figure(1)
imshow(compressed_image,[])
title('Compressed Image')
figure(2)
imshow(hostimage,[])
title('Host Image')

end

% Function to calculate psnr used above
function [p mse]= psnr(original,reconstructed)

Mc= size(original,1);
Nc= size(original,2);
a1= original(1:Mc,1:Nc);
a2= reconstructed(1:Mc,1:Nc);

error=double(a1)-double(a2);

sum=0;
for (i=1:Mc)
for(j=1:Nc)
sum = sum + ((error(i,j))^2);
end
end

mse = sum/(Mc*Nc);
p = 10*log(double((255*255)/mse));

end

5. Digital Watermarking Using Discrete Wavelet Transform (DWT) & Singular Value Decomposition (SVD) & Checking its Robustness in the presence of JPEG Image Compression

function []= DWTSVD_JPEG()

a=10;
% determine size of host image
Mc=size(hostimage,1); %Height
Nc=size(hostimage,2); %Width

% determine size of watermark image
watermark= im2bw(watermark,0.5);

Mm=size(watermark,1); %Height
Nm=size(watermark,2); %Width

[cA1,cH1,cV1,cD1] = dwt2(hostimage,'db1');
figure (1)
imshow(hostimage,[]);
title('Host Image');
figure (2)
subplot(2,2,1);
imshow(cA1,[]);
title('LL');
subplot(2,2,2);
imshow(cH1,[]);
title('LH');
subplot(2,2,3);
imshow(cV1,[]);
title('HL');
subplot(2,2,4);
imshow(cD1,[]);
title('HH');

%Ih= idwt2(cA1,cH1,cV1,cD1,'db1');
Ih= idwt2([],[],[],cD1,'db1');

[Uh,Sh,Vh]= svd(Ih);

[Uw,Sw,Vw]= svd(double(watermark));

Shw= Sh+a*Sw;
VhT=transpose(Vh);
Ihw= Uh*Shw*VhT;
[cA2,cH2,cV2,cD2] = dwt2(Ihw,'db1');
watermarked_image= idwt2(cA1,cH1,cV1,cD2,'db1');

figure(3)
imshow(watermarked_image,[])
title('Watermarked Image');

%//////////////////// JPEG Image Compression ///////////////////

watermarked_image_noisy = jpeg(watermarked_image);
figure(4)
imshow(watermarked_image_noisy,[])
title('Compressed Image')

%//////////////////// Watermark Extraction ////////////////////

[cAm1,cHm1,cVm1,cDm1] = dwt2(watermarked_image_noisy,'db1');

%Ih= idwt2(cA1,cH1,cV1,cD1,'db1');
Ihm= idwt2([],[],[],cDm1,'db1');

[Uhm,Shm,Vhm]= svd(Ihm);

[Uwm,Swm,Vwm]= svd(double(watermark));

Sm= abs((Shm-Swm))/a;

VwmT=transpose(Vwm);

watermark_extracted= Uwm*Sm*VwmT;
watermark_extracted= im2bw(watermark_extracted,0.005);

% display psnr of watermarked image
[p mse] = psnr(hostimage,watermarked_image_noisy)  % Function written below

% display normalized correlation
N= NC(watermark,watermark_extracted)   % Function written below

figure (5)
imshow(watermark,[])
title('Watermark Image');
figure(6)
imshow(watermark_extracted,[])
title('Extracted Watermark Image');

end

% Function to calculate psnr used above
function [p mse]= psnr(original,reconstructed)

Mc= size(original,1);
Nc= size(original,2);
a1= original(1:Mc,1:Nc);
a2= reconstructed(1:Mc,1:Nc);

error=double(a1)-double(a2);

sum=0;
for (i=1:Mc)
for(j=1:Nc)
sum = sum + ((error(i,j))^2);
end
end

mse = sum/(Mc*Nc);
p = 10*log(double((255*255)/mse));

end

% Function used to calculate Normalized Coefficient in above code

function [N]= NC(original,reconstructed)

Mc= size(original,1);
Nc= size(original,2);
sum1 = 0;
sum2 = 0;
sum3 = 0;

for (i=1:Mc)
for(j=1:Nc)

sum1 = sum1 + (original(i,j)*reconstructed(i,j));
sum2 = sum2 + (original(i,j)*original(i,j));
sum3 = sum3 + (reconstructed(i,j)*reconstructed(i,j));

end
end

sqsum2= sqrt(sum2);
sqsum3= sqrt(sum3);

N = (sum1)/(sqsum2*sqsum3);

end