matlab三幅图像融合并计算其fft2频谱

clear;
Irgb_a = imread('1.jpg');    %读入图片
[m, n, l] = size(Irgb_a);
Igray_a = rgb2gray(Irgb_a);    %转化成灰度图
Igfft_a = fft2(Igray_a);        %对灰度图进行dct变换
I1 = im2double(Irgb_a);      %将数据转换成双精度格式
for i = 1:m
    for j = 1:n
        for k = 1:l
            I(i,j+(k-1)*n) = I1(i,j,k);
        end
    end
end
Irfft = fft2(I);
for i = 1:m
    for j = 1:n
        for k = 1:l
            Irifft_a(i,j,k) = I(i,j+(k-1)*n);
        end
    end
end
%imshow(Igray_a);
%title('原始灰度图');
%figure, imshow(log(abs(Igfft_a)),[]), colormap(jet(64)),colorbar;
%title('空间频谱图');
%figure,imshow(Irgb_a);
%title('原始彩图');
%figure,imshow(Irifft_a);
%title('复原彩图');

Irgb_b = imread('2.jpg');    %读入图片
[m, n, l] = size(Irgb_b);
Igray_b = rgb2gray(Irgb_b);    %转化成灰度图
Igfft_b = fft2(Igray_b);        %对灰度图进行dct变换
I1 = im2double(Irgb_b);      %将数据转换成双精度格式
for i = 1:m
    for j = 1:n
        for k = 1:l
            I(i,j+(k-1)*n) = I1(i,j,k);
        end
    end
end
Irfft = fft2(I);
for i = 1:m
    for j = 1:n
        for k = 1:l
            Irifft_b(i,j,k) = I(i,j+(k-1)*n);
        end
    end
end
%imshow(Igray_b);
%title('原始灰度图');
%figure, imshow(log(abs(Igfft_b)),[]), colormap(jet(64)),colorbar;
%title('空间频谱图');
%figure,imshow(Irgb_b);
%title('原始彩图');
%figure,imshow(Irifft_b);
%title('复原彩图');

Irgb_c = imread('3.jpg');    %读入图片
[m, n, l] = size(Irgb_c);
Igray_c = rgb2gray(Irgb_c);    %转化成灰度图
Igfft_c = fft2(Igray_c);        %对灰度图进行dct变换
I1 = im2double(Irgb_c);      %将数据转换成双精度格式
for i = 1:m
    for j = 1:n
        for k = 1:l
            I(i,j+(k-1)*n) = I1(i,j,k);
        end
    end
end
Irfft = fft2(I);
for i = 1:m
    for j = 1:n
        for k = 1:l
            Irifft_c(i,j,k) = I(i,j+(k-1)*n);
        end
    end
end
%imshow(Igray_c);
%title('原始灰度图');
%figure, imshow(log(abs(Igfft_c)),[]), colormap(jet(64)),colorbar;
%title('空间频谱图');
%figure,imshow(Irgb_c);
%title('原始彩图');
%figure,imshow(Irifft_c);
%title('复原彩图');

%叠加部分程序
a=imread('1.jpg');
%a=rgb2gray(a);
b=imread('2.jpg');
%b=rgb2gray(b);
c=imread('3.jpg');
%c=rgb2gray(c);
T=0.25;   %融合的一个比例值
t=[T^2 2*T*(1-T) (1-T)^2];
result_image=a.*t(1)+b.*t(2)+c.*t(3);
%计算频谱
Irgb_r = result_image;    %读入图片
[m, n, l] = size(Irgb_r);
Igray_r = rgb2gray(Irgb_r);    %转化成灰度图
Igfft_r = fft2(Igray_r);        %对灰度图进行dct变换
I1 = im2double(Irgb_r);      %将数据转换成双精度格式
for i = 1:m
    for j = 1:n
        for k = 1:l
            I(i,j+(k-1)*n) = I1(i,j,k);
        end
    end
end
Irfft = fft2(I);
for i = 1:m
    for j = 1:n
        for k = 1:l
            Irifft_r(i,j,k) = I(i,j+(k-1)*n);
        end
    end
end
%imshow(Igray_r);
%title('原始灰度图');
%figure, imshow(log(abs(Igfft_r)),[]), colormap(jet(64)),colorbar;
%title('空间频谱图');
%figure,imshow(Irgb_r);
%title('原始彩图');
%figure,imshow(Irifft_r);
%title('复原彩图');

figure;
subplot(221);
imshow(Irgb_a);
title('原始彩图a');
subplot(222);
imshow(Irgb_b);
title('原始彩图b');
subplot(223);
imshow(Irgb_c);
title('原始彩图c');
subplot(224);
imshow(Irgb_r);
title('原始彩图r');

figure;
subplot(221);
imshow(log(abs(Igfft_a)),[]), colormap(jet(64)),colorbar;
title('空间频谱图a');
subplot(222);
imshow(log(abs(Igfft_b)),[]), colormap(jet(64)),colorbar;
title('空间频谱图b');
subplot(223);
imshow(log(abs(Igfft_c)),[]), colormap(jet(64)),colorbar;
title('空间频谱图c');
subplot(224);
imshow(log(abs(Igfft_r)),[]), colormap(jet(64)),colorbar;
title('空间频谱图r');

figure;
subplot(221);
imshow(log(abs(Igfft_r-Igfft_a)),[]), colormap(jet(64)),colorbar;
title('r频谱-a频谱');
subplot(222);
imshow(log(abs(Igfft_r-Igfft_b)),[]), colormap(jet(64)),colorbar;
title('r频谱-b频谱');
subplot(223);
imshow(log(abs(Igfft_r-Igfft_c)),[]), colormap(jet(64)),colorbar;
title('r频谱-c频谱');
subplot(224);
imshow(log(abs(Igfft_r-Igfft_r)),[]), colormap(jet(64)),colorbar;
title('r频谱-r频谱');

结果例：