1. Performance Analysis of BPSK in the presence of AWGN Noise and Rayleigh Fading Channel
2. Matlab Code for Wiener Filtering of an Image in Frequency Domain
3. Singular Value Decomposition (SVD) Based Image Compression
function []= SVD_C()
hostimage= imread('lena.bmp');
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
hostimage= imread('lena.bmp');
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
function []= DWTSVD_JPEG()
hostimage= imread('lena.bmp');
a=10;
% determine size of host image
Mc=size(hostimage,1); %Height
Nc=size(hostimage,2); %Width
% determine size of watermark image
watermark = imread('watermark_m.bmp');
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
x_bits= randi([0 1],1,1000000); % Input bits
hq = modem.pskmod('M', 2, 'PhaseOffset', 0, 'SymbolOrder', 'Gray', 'InputType', 'Bit');
%Modulation Object for BPSK
%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
% about others or
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()
hostimage= imread('lena.bmp');
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
hostimage= imread('lena.bmp');
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()
hostimage= imread('lena.bmp');
a=10;
% determine size of host image
Mc=size(hostimage,1); %Height
Nc=size(hostimage,2); %Width
% determine size of watermark image
watermark = imread('watermark_m.bmp');
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