Calculating UACI, NPCR for 2 images

Presently I have a problem in implementing the UACI and PSNR code for 2 images which are respectively Original & Cipher Images. I need the code for calculating it......

What would the Code? I don't have much idea about it.. Do I need to find UACI between 2 images or only 1 image will be sufficient

[sajjad]

Professor, i am using a choatic map (in code it is with name a squence, this sequence is gotten from a matrix S2, that S2 (minima of a chaotic system, i have already solved that chatic system no need to solve it again) ) for image encrotion by using algorithms given in this paper (Image encryption using chaos and block cipher, with DOI No10.5539/cis.v4n1p172). But code is not running corecctly. Please checkt it. I have attached the all files. It is an image encryption scheme based on combination of pixel shuffling and new modified version of simplified AES. Chaotic map is used for shuffling and improving S-AES efficiency through S-box design. Chaos is used to expand diffusion and confusion in the image.

'' in the paper he use a baker's map "

(FOLLOWIGN IS THE CODE FOR MY MAP i,e SEQUENCE instead OF BAKER's MAP)

%% (1)

% Function for Shuffling Pixels Using Chaotic Sequence (Replacing Baker's Map)

function shuffled_img = shuffle_pixels_with_sequence(img, sequence)

[rows, cols] = size(img);

num_points = size(sequence, 1);

% Generate initial grid coordinates

[X, Y] = meshgrid(1:cols, 1:rows);

X = X(:) / cols;

Y = Y(:) / rows;

% Iterate to compute the shuffled positions using the chaotic sequence

for i = 1:num_points

X_next = mod(X + sequence(i, 1), 1);

Y_next = mod(Y + sequence(i, 2), 1);

X = X_next;

Y = Y_next;

end

% Ensure X and Y values are within the valid range

X = ceil(X * cols);

Y = ceil(Y * rows);

% Make sure X and Y are within the range [1, rows] and [1, cols]

X = max(min(X, cols), 1);

Y = max(min(Y, rows), 1);

% Rearrange pixels based on shuffled indices

indices = sub2ind([rows, cols], Y, X);

shuffled_img = img(indices);

shuffled_img = reshape(shuffled_img, rows, cols);

% **Step 1: Vertical Permutation**

Pmap = randperm(rows);

shuffled_img = shuffled_img(Pmap, :);

% **Step 2: Horizontal Permutation**

Pmap = randperm(cols);

shuffled_img = shuffled_img(:, Pmap);

% **Step 3: Apply Circular Shift for Diagonal Permutation**

for it = 2:rows

Shiftsize = mod(num_points - it + 1, cols); % Prevents exceeding bounds

shuffled_img(it, :) = circshift(shuffled_img(it, :), [0, Shiftsize]);

end

% **Step 4: Apply Vertical Permutation Again**

shuffled_img = shuffled_img(Pmap, :);

% **Step 5: Apply Final Diagonal Permutation**

for it = 2:rows

Shiftsize = mod(num_points - it + 1, cols);

shuffled_img(it, :) = circshift(shuffled_img(it, :), [0, -Shiftsize]);

end

end

%% (2)

function encrypted_img = s_aes_encrypt(img, sbox)

[rows, cols] = size(img);

% Ensure the S-box is expanded to match the image size

if size(sbox, 1) < rows || size(sbox, 2) < cols

sbox_expanded = repmat(sbox, ceil(rows / 256), ceil(cols / 256));

sbox_expanded = sbox_expanded(1:rows, 1:cols);

else

sbox_expanded = sbox(1:rows, 1:cols);

end

% Convert image to double for computations

img = double(img);

% Round 1

img = add_round_key(img, sbox_expanded); % Add Round Key

img = substitute_nibbles(img, sbox_expanded); % Substitute Nibbles

img = shift_rows(img); % Shift Rows

img = mix_columns(img); % Mix Columns

img = add_round_key(img, sbox_expanded); % Add Round Key

% Round 2

img = substitute_nibbles(img, sbox_expanded); % Substitute Nibbles

img = shift_rows(img); % Shift Rows

img = add_round_key(img, sbox_expanded); % Add Round Key

% Convert back to uint8 for image display

encrypted_img = uint8(img);

end

%% (3)

function state = substitute_nibbles(state, sbox)

[rows, cols] = size(state);

% Substitute each nibble using the S-box

for i = 1:rows

for j = 1:cols

% Ensure the value is within the range of the S-box

nibble_value = mod(state(i, j), 256) + 1; % Map to 1-256

state(i, j) = sbox(nibble_value);

end

end

end

%% (4)

function state = shift_rows(state)

% Shift rows as per S-AES

state(2, :) = circshift(state(2, :), -1); % Shift second row by 1 position

end

%%% (5)

function state = mix_columns(state)

% Mix columns as per S-AES

temp = state;

state(1, 1) = bitxor(temp(1, 1), temp(2, 1));

state(1, 2) = bitxor(temp(1, 2), temp(2, 2));

state(2, 1) = bitxor(temp(1, 1), temp(2, 1));

state(2, 2) = bitxor(temp(1, 2), temp(2, 2));

end

%%% (6)

function state = add_round_key(state, round_key)

% XOR the state with the round key

state = bitxor(state, round_key);

end

%%%%(7)

%% (5) generate_sbox - Fully Corrected Version

function sbox = generate_sbox(sequence)

NoIt = length(sequence);

D = zeros(1, NoIt);

% **Step 1-5: Compute D(it) Using Iterative Transformation**

for it = 1:NoIt

idx = mod(it - 1, length(sequence)) + 1;

X = sequence(idx, 1);

Y = sequence(idx, 2);

% ? Corrected: Use iterative transformations

X_next = mod(X + Y, 1);

Y_next = mod(X + 2 * Y, 1);

% Update X, Y for next iteration

X = X_next;

Y = Y_next;

% ? Compute D(it) using correct transformation

D(it) = sqrt(abs(X.^3 + Y.^3));

end

% **Step 6-10: Apply Sorting & Nonlinear Transformation**

D = sort(D, 'ascend');

M = max(D);

for it1 = 1:NoIt

for it2 = 1:NoIt

if D(it2) == M

D(it2) = -16 + it1;

end

end

end

% **Step 11-20: Normalize & Create the S-Box**

D = abs(D);

max_D = max(D);

sbox = uint8(mod(round(D * 255 / max_D), 256));

% **Ensure the S-Box is exactly (256 × 256)**

sbox = repmat(sbox(:), ceil(256*256/numel(sbox)), 1);

sbox = reshape(sbox(1:256*256), 256, 256);

end

%%%%(8)

%% (8) encryption_main_subscript - Fully Matches Figure 5

clc; clear; close all;

%% **Step 1: Load and Preprocess Image**

filePath = 'C:\Users\HP EliteBook 840 G3\Desktop\pic3.jpeg';

img = imread(filePath);

% Convert to grayscale if the image is RGB

if size(img, 3) == 3

img = rgb2gray(img);

end

% Resize image to 256x256 (if not already)

img = imresize(img, [256, 256]);

% Display (a) Original Image

figure;

subplot(2, 3, 1);

imshow(img);

title('(a) Original Image');

%% **Step 2: Load Sequence from S2.xlsx for Pixel Shuffling**

sequenceFile = 'C:\Users\HP EliteBook 840 G3\Desktop\S1.xlsx';

data = xlsread(sequenceFile);

data = data(:); % Flatten matrix

% Scale the sequence between 0 and 1 (Normalization)

scaled_sequence = (data - min(data)) / (max(data) - min(data));

maxima_x = scaled_sequence(1:end-1);

maxima_y = scaled_sequence(2:end);

sequence = [maxima_x, maxima_y];

%% **Step 3: Apply Algorithm 1 (Pixel Shuffling)**

shuffledImg = shuffle_pixels_with_sequence(img, sequence);

% Display (b) Shuffled Image

subplot(2, 3, 2);

imshow(shuffledImg);

title('(b) Shuffled Image');

%% **Step 4: Apply Algorithm 2 (Generate Key-Dependent S-Box)**

sbox = generate_sbox(sequence); % Generate dynamic S-Box

%% **Step 5: Apply Simplified AES (S-AES) Encryption**

I = s_aes_encrypt(shuffledImg, sbox); % Pass S-Box explicitly

% Display (c) Ciphered Image

subplot(2, 3, 3);

imshow(I);

title('(c) Encrypted Image');

%% **Step 6: Compute and Display Histograms**

subplot(2, 3, 4); imhist(img); title('(d) Histogram of Original Image');

subplot(2, 3, 5); imhist(shuffledImg); title('(e) Histogram of Shuffled Image');

subplot(2, 3, 6); imhist(I); title('(f) Histogram of Encrypted Image');

% %% **Step 7: Save Encrypted Image**

% outputPath = 'C:\Users\HP EliteBook 840 G3\Desktop\encrypted_image.png';

% imwrite(cipheredImg, outputPath);

% disp(['Encrypted image saved to: ', outputPath]);

[sajjad]

For S1; from the following code just consider MAXIMA and in above codes i have converted S1 in sequence.

The system this paper (DOI: 10.1017/S0022112003003835) Thank you.

%% funn

function dydt = funn(t,y,x,nx,dx,delta,Reynolds,H,Hdot);

%persistent iter

%if isempty(iter)

% iter = 0;

%end

%iter = iter + 1;

%if mod(iter,1000)==0

% t

% iter = 0;

%end

G = y(1:nx);

F = y(nx+1:2*nx);

% F

% Compute spatial derivatives

dFdx = zeros(nx,1);

d2Fdx2 = zeros(nx,1);

dFdx(2:nx-1) = (F(3:nx)-F(1:nx-2))/(2*dx);

dFdx(nx) = Hdot(t); % This corresponds to F'(1,t) = H'(t)

d2Fdx2(2:nx-1) = (F(3:nx)-2*F(2:nx-1)+F(1:nx-2))/dx^2;

%d2Fdx2(nx) = (dFdx(nx)-dFdx(nx-1))/dx;

Fnxp1 = F(nx-1) + 2*dx*dFdx(nx);

d2Fdx2(nx) = (Fnxp1-2*F(nx)+F(nx-1))/dx^2;

% Compute temporal derivatives

dFdt = zeros(nx,1);

dFdt(1) = F(1);

dFdt(2:nx-1) = d2Fdx2(2:nx-1)+dFdx(2:nx-1)./x(2:nx-1)-...

F(2:nx-1)./x(2:nx-1).^2+H(t)^2*G(2:nx-1);

dFdt(nx) = F(nx) + Hdot(t); % This corresponds to F(1,t) =- H'(t)

% G

% Compute spatial derivatives

dGdx = zeros(nx,1);

d2Gdx2 = zeros(nx,1);

dGdx(2:nx-1) = (G(3:nx)-G(1:nx-2))/(2*dx);

d2Gdx2(2:nx-1) = (G(3:nx)-2*G(2:nx-1)+G(1:nx-2))/dx^2;

% Compute temporal derivatives

dGdt = zeros(nx,1);

dGdt(1) = G(1);

dGdt(2:nx-1) = Hdot(t)/H(t).*x(2:nx-1).*dGdx(2:nx-1)+ ...

1/H(t).*F(2:nx-1).*dGdx(2:nx-1)- ...

1/H(t).*dFdx(2:nx-1).*G(2:nx-1)- ...

2.*F(2:nx-1).*G(2:nx-1)./(H(t)*x(2:nx-1))+ ...

(d2Gdx2(2:nx-1)+dGdx(2:nx-1)./x(2:nx-1)-G(2:nx-1)./x(2:nx-1).^2)./ ...

(H(t)^2*Reynolds);

dGdt(nx) = G(nx) - (d2Fdx2(nx)+dFdx(nx)/x(nx)-F(nx)/x(nx)^2)/(-H(t)^2);

%Taken from the article, should be the same as set

%dGdt(nx) = G(nx) + 2/H(t)^2*((F(nx-1)+Hdot(t)*(1+dx))/dx^2 + Hdot(t));

dydt = [dGdt;dFdt];

end

%%% MAIN SUBSCRIPTT

format long;

clear all;

close all;

clc;

% Parameters

delta = 0.3; % Amplitude of H

Reynolds = 900; % Reynolds number

tstart = 0; % Start time

tend = 206400; % End time

nx = 101; % Number of spatial points

xstart = 0.0; % Spatial domain start

xend = 1.0; % Spatial domain end

x = linspace(xstart, xend, nx).'; % Spatial grid

dx = (xend - xstart) / (nx - 1); % Spatial step size

tspan = [tstart, tend]; % Time span

% numpoints = 1000;

% tspan = linspace(tstart, tend,numpoints); % Time span

% Initial conditions

G0 = zeros(nx, 1); % Initial condition for G

F0 = zeros(nx, 1); % Initial condition for F

y0 = [G0; F0]; % Combine into one vector

% Define H(t) and its derivative

H = @(t) 1 + delta * cos(2 * t);

Hdot = @(t) -2 * delta * sin(2 * t);

% Mass matrix for DAE

M = [eye(nx), zeros(nx); zeros(nx, 2*nx)];

M(1, 1) = 0; % Boundary condition G(0,t) = 0

M(nx, nx) = 0; % Boundary condition G(1,t) = 0

options = odeset('Mass', M, 'RelTol', 1e-3, 'AbsTol', 1e-6);

% Solve the system

[T, Y] = ode23t(@(t, y) funn(t, y, x, nx, dx, delta, Reynolds, H, Hdot), tspan, y0, options);

% Extract G and F

G = Y(:, 1:nx);

F = Y(:, nx+1:2*nx);

% Evaluate G at eta = 1 using interpolation

eta_query = 1; % Change from 1/4 to 1

G_eta_1 = zeros(size(T)); % Preallocate

for i = 1:length(T)

G_eta_1(i) = interp1(x, G(i, :), eta_query, 'spline'); % Interpolation

end

% Plot G(1, t) over time

figure(1);

plot(T, G_eta_1, 'LineWidth', 1.5);

xlabel('Time (t)');

ylabel('G(1, t)');

title('Time evolution of G(1, t)');

grid on;

[maxima_G_eta_1, locs_max_G_eta_1] = findpeaks(G_eta_1); % Maxima of G(1/4, t)

[minima_G_eta_1, locs_min_G_eta_1] = findpeaks(-G_eta_1); % Minima of G(1/4, t)

minima_G_eta_1 = -minima_G_eta_1; % Correct sign of minima

% Display counts of maxima and minima

sajjad1 = length(maxima_G_eta_1);

sajjad2 = length(minima_G_eta_1);

disp(['Number of maxima: ', num2str(sajjad1)]);

disp(['Number of minima: ', num2str(sajjad2)]);

% Multiply maxima and minima by 10^15 and apply mod 256

scaled_maxima = mod(round(maxima_G_eta_1 * 1e15), 256);

scaled_minima = mod(round(minima_G_eta_1 * 1e15), 256);

% Display scaled maxima and minima

disp('Scaled maxima (mod 256):');

disp(scaled_maxima);

disp('Scaled minima (mod 256):');

disp(scaled_minima);

% Return Map for Maxima

figure(2); % New figure

plot(maxima_G_eta_1(1:end-1), maxima_G_eta_1(2:end), 'o','LineWidth', 2,'Color', 'r', 'MarkerSize', 4);

xlabel('Maxima_{n}');

ylabel('Maxima_{n+1}');

title('Return Map of Maxima');

grid on;

% Return Map for Minima

figure(3); % Another new figure

plot(minima_G_eta_1(1:end-1), minima_G_eta_1(2:end), 'o','LineWidth', 2,'Color', 'g', 'MarkerSize', 4);

xlabel('Minima_{n}');

ylabel('Minima_{n+1}');

title('Return Map of Minima');

grid on;

Please feel to ask any other information about it.

Thank you.

[sajjad]

Pic1 one is for encryption.

[sajjad]

Professor, i am using a choatic map (in code it is with name a squence, this sequence is gotten from a matrix S2, that S2 (minima of a chaotic system, i have already solved that chatic system no need to solve it again) ) for image encrotion by using algorithms given in this paper (Image encryption using chaos and block cipher, with DOI No10.5539/cis.v4n1p172). But code is not running corecctly. Please checkt it. I have attached the all files. It is an image encryption scheme based on combination of pixel shuffling and new modified version of simplified AES. Chaotic map is used for shuffling and improving S-AES efficiency through S-box design. Chaos is used to expand diffusion and confusion in the image.

'' in the paper he use a baker's map "

(FOLLOWIGN IS THE CODE FOR MY MAP i,e SEQUENCE instead OF BAKER's MAP)

%% (1)

% Function for Shuffling Pixels Using Chaotic Sequence (Replacing Baker's Map by my sequence)

function shuffled_img = shuffle_pixels_with_sequence(img, sequence)

[rows, cols] = size(img);

num_points = size(sequence, 1);

% Generate initial grid coordinates

[X, Y] = meshgrid(1:cols, 1:rows);

X = X(:) / cols;

Y = Y(:) / rows;

% Iterate to compute the shuffled positions using the chaotic sequence

for i = 1:num_points

X_next = mod(X + sequence(i, 1), 1);

Y_next = mod(Y + sequence(i, 2), 1);

X = X_next;

Y = Y_next;

end

% Ensure X and Y values are within the valid range

X = ceil(X * cols);

Y = ceil(Y * rows);

% Make sure X and Y are within the range [1, rows] and [1, cols]

X = max(min(X, cols), 1);

Y = max(min(Y, rows), 1);

% Rearrange pixels based on shuffled indices

indices = sub2ind([rows, cols], Y, X);

shuffled_img = img(indices);

shuffled_img = reshape(shuffled_img, rows, cols);

% **Step 1: Vertical Permutation**

Pmap = randperm(rows);

shuffled_img = shuffled_img(Pmap, :);

% **Step 2: Horizontal Permutation**

Pmap = randperm(cols);

shuffled_img = shuffled_img(:, Pmap);

% **Step 3: Apply Circular Shift for Diagonal Permutation**

for it = 2:rows

Shiftsize = mod(num_points - it + 1, cols); % Prevents exceeding bounds

shuffled_img(it, :) = circshift(shuffled_img(it, :), [0, Shiftsize]);

end

% **Step 4: Apply Vertical Permutation Again**

shuffled_img = shuffled_img(Pmap, :);

% **Step 5: Apply Final Diagonal Permutation**

for it = 2:rows

Shiftsize = mod(num_points - it + 1, cols);

shuffled_img(it, :) = circshift(shuffled_img(it, :), [0, -Shiftsize]);

end

end

%% (2)

function encrypted_img = s_aes_encrypt(img, sbox)

[rows, cols] = size(img);

% Ensure the S-box is expanded to match the image size

if size(sbox, 1) < rows || size(sbox, 2) < cols

sbox_expanded = repmat(sbox, ceil(rows / 256), ceil(cols / 256));

sbox_expanded = sbox_expanded(1:rows, 1:cols);

else

sbox_expanded = sbox(1:rows, 1:cols);

end

% Convert image to double for computations

img = double(img);

% Round 1

img = add_round_key(img, sbox_expanded); % Add Round Key

img = substitute_nibbles(img, sbox_expanded); % Substitute Nibbles

img = shift_rows(img); % Shift Rows

img = mix_columns(img); % Mix Columns

img = add_round_key(img, sbox_expanded); % Add Round Key

% Round 2

img = substitute_nibbles(img, sbox_expanded); % Substitute Nibbles

img = shift_rows(img); % Shift Rows

img = add_round_key(img, sbox_expanded); % Add Round Key

% Convert back to uint8 for image display

encrypted_img = uint8(img);

end

%% (3)

function state = substitute_nibbles(state, sbox)

[rows, cols] = size(state);

% Substitute each nibble using the S-box

for i = 1:rows

for j = 1:cols

% Ensure the value is within the range of the S-box

nibble_value = mod(state(i, j), 256) + 1; % Map to 1-256

state(i, j) = sbox(nibble_value);

end

end

end

%% (4)

function state = shift_rows(state)

% Shift rows as per S-AES

state(2, :) = circshift(state(2, :), -1); % Shift second row by 1 position

end

%% (5)

function state = mix_columns(state)

% Mix columns as per S-AES

temp = state;

state(1, 1) = bitxor(temp(1, 1), temp(2, 1));

state(1, 2) = bitxor(temp(1, 2), temp(2, 2));

state(2, 1) = bitxor(temp(1, 1), temp(2, 1));

state(2, 2) = bitxor(temp(1, 2), temp(2, 2));

end

%% (6)

function state = add_round_key(state, round_key)

% XOR the state with the round key

state = bitxor(state, round_key);

end

%% (7)

%% (5) generate_sbox - Fully Corrected Version

function sbox = generate_sbox(sequence)

NoIt = length(sequence);

D = zeros(1, NoIt);

% **Step 1-5: Compute D(it) Using Iterative Transformation**

for it = 1:NoIt

idx = mod(it - 1, length(sequence)) + 1;

X = sequence(idx, 1);

Y = sequence(idx, 2);

% ? Corrected: Use iterative transformations

X_next = mod(X + Y, 1);

Y_next = mod(X + 2 * Y, 1);

% Update X, Y for next iteration

X = X_next;

Y = Y_next;

% ? Compute D(it) using correct transformation

D(it) = sqrt(abs(X.^3 + Y.^3));

end

% **Step 6-10: Apply Sorting & Nonlinear Transformation**

D = sort(D, 'ascend');

M = max(D);

for it1 = 1:NoIt

for it2 = 1:NoIt

if D(it2) == M

D(it2) = -16 + it1;

end

end

end

% **Step 11-20: Normalize & Create the S-Box**

D = abs(D);

max_D = max(D);

sbox = uint8(mod(round(D * 255 / max_D), 256));

% **Ensure the S-Box is exactly (256 × 256)**

sbox = repmat(sbox(:), ceil(256*256/numel(sbox)), 1);

sbox = reshape(sbox(1:256*256), 256, 256);

end

%%% (8)

%% (8) encryption_main_subscript - Fully Matches Figure 5

clc; clear; close all;

%% **Step 1: Load and Preprocess Image**

filePath = 'C:\Users\HP EliteBook 840 G3\Desktop\pic1.jpeg';

img = imread(filePath);

% Convert to grayscale if the image is RGB

if size(img, 3) == 3

img = rgb2gray(img);

end

% Resize image to 256x256 (if not already)

img = imresize(img, [256, 256]);

% Display (a) Original Image

figure;

subplot(2, 3, 1);

imshow(img);

title('(a) Original Image');

%% **Step 2: Load Sequence from S2.xlsx for Pixel Shuffling**

sequenceFile = 'C:\Users\HP EliteBook 840 G3\Desktop\S2.xlsx';

data = xlsread(sequenceFile);

data = data(:); % Flatten matrix

% Scale the sequence between 0 and 1 (Normalization)

scaled_sequence = (data - min(data)) / (max(data) - min(data));

maxima_x = scaled_sequence(1:end-1);

maxima_y = scaled_sequence(2:end);

sequence = [maxima_x, maxima_y];

%% **Step 3: Apply Algorithm 1 (Pixel Shuffling)**

shuffledImg = shuffle_pixels_with_sequence(img, sequence);

% Display (b) Shuffled Image

subplot(2, 3, 2);

imshow(shuffledImg);

title('(b) Shuffled Image');

%% **Step 4: Apply Algorithm 2 (Generate Key-Dependent S-Box)**

sbox = generate_sbox(sequence); % Generate dynamic S-Box

%% **Step 5: Apply Simplified AES (S-AES) Encryption**

I = s_aes_encrypt(shuffledImg, sbox); % Pass S-Box explicitly

% Display (c) Ciphered Image

subplot(2, 3, 3);

imshow(I);

title('(c) Encrypted Image');

%% **Step 6: Compute and Display Histograms**

subplot(2, 3, 4); imhist(img); title('(d) Histogram of Original Image');

subplot(2, 3, 5); imhist(shuffledImg); title('(e) Histogram of Shuffled Image');

subplot(2, 3, 6); imhist(I); title('(f) Histogram of Encrypted Image');

% %% **Step 7: Save Encrypted Image**

% outputPath = 'C:\Users\HP EliteBook 840 G3\Desktop\encrypted_image.png';

% imwrite(cipheredImg, outputPath);

% disp(['Encrypted image saved to: ', outputPath]);

[sajjad]

For S1; from the following code just consider MAXIMA and in above codes i have converted S1 in sequence.

The system this paper (DOI: 10.1017/S0022112003003835) Thank you.

%% funn

function dydt = funn(t,y,x,nx,dx,delta,Reynolds,H,Hdot);

%persistent iter

%if isempty(iter)

% iter = 0;

%end

%iter = iter + 1;

%if mod(iter,1000)==0

% t

% iter = 0;

%end

G = y(1:nx);

F = y(nx+1:2*nx);

% F

% Compute spatial derivatives

dFdx = zeros(nx,1);

d2Fdx2 = zeros(nx,1);

dFdx(2:nx-1) = (F(3:nx)-F(1:nx-2))/(2*dx);

dFdx(nx) = Hdot(t); % This corresponds to F'(1,t) = H'(t)

d2Fdx2(2:nx-1) = (F(3:nx)-2*F(2:nx-1)+F(1:nx-2))/dx^2;

%d2Fdx2(nx) = (dFdx(nx)-dFdx(nx-1))/dx;

Fnxp1 = F(nx-1) + 2*dx*dFdx(nx);

d2Fdx2(nx) = (Fnxp1-2*F(nx)+F(nx-1))/dx^2;

% Compute temporal derivatives

dFdt = zeros(nx,1);

dFdt(1) = F(1);

dFdt(2:nx-1) = d2Fdx2(2:nx-1)+dFdx(2:nx-1)./x(2:nx-1)-...

F(2:nx-1)./x(2:nx-1).^2+H(t)^2*G(2:nx-1);

dFdt(nx) = F(nx) + Hdot(t); % This corresponds to F(1,t) =- H'(t)

% G

% Compute spatial derivatives

dGdx = zeros(nx,1);

d2Gdx2 = zeros(nx,1);

dGdx(2:nx-1) = (G(3:nx)-G(1:nx-2))/(2*dx);

d2Gdx2(2:nx-1) = (G(3:nx)-2*G(2:nx-1)+G(1:nx-2))/dx^2;

% Compute temporal derivatives

dGdt = zeros(nx,1);

dGdt(1) = G(1);

dGdt(2:nx-1) = Hdot(t)/H(t).*x(2:nx-1).*dGdx(2:nx-1)+ ...

1/H(t).*F(2:nx-1).*dGdx(2:nx-1)- ...

1/H(t).*dFdx(2:nx-1).*G(2:nx-1)- ...

2.*F(2:nx-1).*G(2:nx-1)./(H(t)*x(2:nx-1))+ ...

(d2Gdx2(2:nx-1)+dGdx(2:nx-1)./x(2:nx-1)-G(2:nx-1)./x(2:nx-1).^2)./ ...

(H(t)^2*Reynolds);

dGdt(nx) = G(nx) - (d2Fdx2(nx)+dFdx(nx)/x(nx)-F(nx)/x(nx)^2)/(-H(t)^2);

%Taken from the article, should be the same as set

%dGdt(nx) = G(nx) + 2/H(t)^2*((F(nx-1)+Hdot(t)*(1+dx))/dx^2 + Hdot(t));

dydt = [dGdt;dFdt];

end

%%% MAIN SUBSCRIPTT

format long;

clear all;

close all;

clc;

% Parameters

delta = 0.3; % Amplitude of H

Reynolds = 900; % Reynolds number

tstart = 0; % Start time

tend = 206400; % End time

nx = 101; % Number of spatial points

xstart = 0.0; % Spatial domain start

xend = 1.0; % Spatial domain end

x = linspace(xstart, xend, nx).'; % Spatial grid

dx = (xend - xstart) / (nx - 1); % Spatial step size

tspan = [tstart, tend]; % Time span

% numpoints = 1000;

% tspan = linspace(tstart, tend,numpoints); % Time span

% Initial conditions

G0 = zeros(nx, 1); % Initial condition for G

F0 = zeros(nx, 1); % Initial condition for F

y0 = [G0; F0]; % Combine into one vector

% Define H(t) and its derivative

H = @(t) 1 + delta * cos(2 * t);

Hdot = @(t) -2 * delta * sin(2 * t);

% Mass matrix for DAE

M = [eye(nx), zeros(nx); zeros(nx, 2*nx)];

M(1, 1) = 0; % Boundary condition G(0,t) = 0

M(nx, nx) = 0; % Boundary condition G(1,t) = 0

options = odeset('Mass', M, 'RelTol', 1e-3, 'AbsTol', 1e-6);

% Solve the system

[T, Y] = ode23t(@(t, y) funn(t, y, x, nx, dx, delta, Reynolds, H, Hdot), tspan, y0, options);

% Extract G and F

G = Y(:, 1:nx);

F = Y(:, nx+1:2*nx);

% Evaluate G at eta = 1 using interpolation

eta_query = 1; % Change from 1/4 to 1

G_eta_1 = zeros(size(T)); % Preallocate

for i = 1:length(T)

G_eta_1(i) = interp1(x, G(i, :), eta_query, 'spline'); % Interpolation

end

% Plot G(1, t) over time

figure(1);

plot(T, G_eta_1, 'LineWidth', 1.5);

xlabel('Time (t)');

ylabel('G(1, t)');

title('Time evolution of G(1, t)');

grid on;

[maxima_G_eta_1, locs_max_G_eta_1] = findpeaks(G_eta_1); % Maxima of G(1/4, t)

[minima_G_eta_1, locs_min_G_eta_1] = findpeaks(-G_eta_1); % Minima of G(1/4, t)

minima_G_eta_1 = -minima_G_eta_1; % Correct sign of minima

% Display counts of maxima and minima

sajjad1 = length(maxima_G_eta_1);

sajjad2 = length(minima_G_eta_1);

disp(['Number of maxima: ', num2str(sajjad1)]);

disp(['Number of minima: ', num2str(sajjad2)]);

% Multiply maxima and minima by 10^15 and apply mod 256

scaled_maxima = mod(round(maxima_G_eta_1 * 1e15), 256);

scaled_minima = mod(round(minima_G_eta_1 * 1e15), 256);

% Display scaled maxima and minima

disp('Scaled maxima (mod 256):');

disp(scaled_maxima);

disp('Scaled minima (mod 256):');

disp(scaled_minima);

% Return Map for Maxima

figure(2); % New figure

plot(maxima_G_eta_1(1:end-1), maxima_G_eta_1(2:end), 'o','LineWidth', 2,'Color', 'r', 'MarkerSize', 4);

xlabel('Maxima_{n}');

ylabel('Maxima_{n+1}');

title('Return Map of Maxima');

grid on;

% Return Map for Minima

figure(3); % Another new figure

plot(minima_G_eta_1(1:end-1), minima_G_eta_1(2:end), 'o','LineWidth', 2,'Color', 'g', 'MarkerSize', 4);

xlabel('Minima_{n}');

ylabel('Minima_{n+1}');

title('Return Map of Minima');

grid on;

Please feel to ask any other information about it.

Thank you.

Answer 1

[Image Analyst]

See my demo:

% Demo to calculate PSNR of a gray scale image.
% http://en.wikipedia.org/wiki/PSNR
% by ImageAnalyst 
clc; 
close all; 
clearvars; 
workspace;
% Read in standard MATLAB demo image.
grayImage = imread('cameraman.tif');
[rows columns] = size(grayImage);
subplot(2, 2, 1);
imshow(grayImage, []);
title('Original Grey Scale Image');
set(gcf, 'Position', get(0,'Screensize')); % Maximize figure.
% Add noise to it.
noisyImage = imnoise(grayImage, 'gaussian', 0, 0.003);
subplot(2, 2, 2);
imshow(noisyImage, []);
title('Noisy Image');

% Calculate mean square error.
mseImage = (double(grayImage) - double(noisyImage)) .^ 2;
subplot(2, 2, 3);
imshow(mseImage, []);
title('MSE Image');
mse = sum(sum(mseImage)) / (rows * columns);
% Calculate PSNR (Peak Signal to noise ratio).
PSNR = 10 * log10( 256^2 / mse);
message = sprintf('The mean square error is %.2f\nThe PSNR = %.2f', ...
  mse, PSNR);
msgbox(message);
% set(gcf, 'Position', get(0,'Screensize')); % Maximize figure.

After this example you should be able to use similar code to easily calculate NPCR and UACI, according to http://www.cyberjournals.com/Papers/Apr2011/05.pdf. You need both images to calculate them since they are a measure of the differences between the two images.

[Komeil]

The NPCR is different from the PSNR.

[Walter Roberson]

Image Analyst says "After this example you should be able to use similar code to easily calculate NPCR" -- in other words, he has not given the code here, but he has given a framework that can be adapted to calculate NPCR by anyone who looks up the formula for NPCR.

Answer 2

[Hyderkkk]

How do i get the NPCR and UACI of an original speech and encryption in matlab code ?

[Image Analyst]

Check the File Exchange: NPCR and UACI measurements with statistical tests

[Hyderkkk]

This file give the relation between two encrypted images of same size and type but I need the relation between original speech and encrypted speech . thank you for answer

[Image Analyst]

You're welcome, but that's all I can contribute. If that doesn't work for you, you'll have to start writing code.

[merr LAMINe]

if you use MATLAB , the samples of the speech normalized in the interval [-1,1], so as a personal idea, I think it calculated with this code (Rmrq : uaci is calculated between 2 encrypted speech )

function [ uc_value1 ] = uaci( s1,s2 )

l=length(s1);
sm=sum(abs(s1-s2));
uc_value1=(100/l)*(sm/2);%2 instead of 65535 (1-(-1)=2)

end

[Kiki Kiosk]

Hey, can you help me?

I must calculate uaci between two audio files, but it's not working right ..

This is my code

uaci_score_mine = sum(abs(data(:) - data1(:))) / 2/ N*100;

I have tried with

uaci_score_mine = normcdf(sum(abs(data(:) - data1(:))) / 2/ N)*100;

and

uaci_score_mine = normcdf(sum(abs(data(:) - data1(:))) / 65535/ N)*100;

and

uaci_score_mine = normcdf(sum(abs(data(:) - data1(:)) / 65535)/ N)*100;

but the results are still not that i expect ..

data and data1 are the audio files

N - length(data);

10x for your attention

[Walter Roberson]

Not clear to me where the 65535 is coming from? If you had used 65536 I might have thought it was the same as the 256^2 PSNR formula that Image Analyst used ?

[Kiki Kiosk]

This is the formula I used to calculate.

[Walter Roberson]

are you reading the audio with the 'native' option?

[Kiki Kiosk]

this is my code:

[data,fs]=audioread('Sample1.wav','native');
[data1,fs2]=audioread('SampleCrypt.wav','native');
 N = length(data);
uaci_score_mine = normcdf(sum(abs(data(:) - data1(:))) / 2/ N)*100; 
fprintf('\n The uaci value is %0.15f \n', uaci_score_mine);

[Samaa Yasser]

@Image Analyst hey, please i want to ask question, i want to know why should the value of uaci in image encryption is about 33. ?

[Image Analyst]

I have no idea. Who said uaci should always be about 33?

[Walter Roberson]

This file give the relation between two encrypted images of same size and type but I need the relation between original speech and encrypted speech

Are you familiar with the encryption technique that is often called ROT13? It is a caeser cypher in which each letter in the first half of the alphabet is replaced with the letter in the corresponding position in the second half of the alphabet, and vice versa, so 'A' gets replaced with the character 13 further on (which is 'N'), and 'p' would get replaced with the character 13 earlier (which is 'c')

Now take the original image, and apply ROT13 to it. Then take the result, and apply ROT13 to that, for even stronger encryption.

Now take the relationship between that double-ROT13 encrypted image and the encryption you found using a different method. As they are both encrypted image, there should be no problem comparing the two using an algorithm that compares encrypted images, right?

... But it turns out that applying ROT13 twice always gets you back exactly the original data. Therefore the original data is itself an encrypted version of the original data... just encrypted with a very weak encryption. And so you can use your algorithm that compares two encrypted images in order to compare the original image to the version of it you encrypted some other way than double-ROT13.

[sajjad]

I have almost completed the image encryption. Now how do i calculate the NPCR and UACI of an original image and encryption in matlab code ?

Please guide me in this regard.

Thank you.

[Image Analyst]

@sajjad, did you see my comment above: https://www.mathworks.com/matlabcentral/answers/36171-calculating-uaci-npcr-for-2-images#comment_991829

That's about all I can offer, is that File Exchange entry. I don't have any code for it myself.