FFT求频谱图和功率谱密度图

频谱图

声音频率与能量的关系用频谱表示。在实际使用中，频谱图有三种，即线性振幅谱、对数振幅谱、自功率谱。

线性振幅谱的纵坐标有明确的物理量纲，是最常用的。

对数振幅谱中各谱线的振幅都作了对数计算，所以其纵坐标的单位是dB（分贝）。这个变换的目的是使那些振幅较低的成分相对高振幅成分得以拉高，以便观察掩盖在低幅噪声中的周期信号。

自功率谱是先对测量信号作自相关卷积，目的是去掉随机干扰噪声，保留并突出周期性信号，损失了相位特征，然后再作傅里叶变换。自功率谱图使得周期性信号更加突出。

功率谱图

功率谱是功率谱密度函数的简称，它定义为单位频带内的信号功率。它表示了信号功率随着频率的变化情况，即信号功率在频域的分布状况。

功率谱表示了信号功率随着频率的变化关系。常用于功率信号（区别于能量信号）的表述与分析，其曲线（即功率谱曲线）一般横坐标为频率，纵坐标为功率。

由于功率没有负值，所以功率谱曲线上的纵坐标也没有负数值，功率谱曲线所覆盖的面积在数值上等于信号的总功率（能量）。

时域和频域能量相等。

Parseval定理



有限上序列x{k}的离散fourier变换是正交变换，满足Parseval能量守恒定理，反映了序列在时域的能量等于其变换域的能量。

关于能量定义：信号幅度平方的积分，如果是数字信号，能量就是各点信号幅度值平方后的求和。

论坛帖子中关于等式关系给出的结论是：

求和 (x(tn)^2)T=RMS^2*Ttotal=求和(P(fn))△f*Ttotal

其中，x(tn)是n个x(t)时域采样数据，T是时间间隔，Ttotal是时间总长，P(fn)是第n个功率谱密度值，△f是FFT频率间隔.

最后的结论是相等的，但是信号的能量到底是sum(x.^2)，还是sum(x.^2)*T？按照定义来说是前者没错。但是绝对的能量计算若不跟采样频率（采样间隔）结合起来，又有什么对比作用？

同样1000个点幅值为1，一组波形是1秒内采到的，另一组波形是10秒内采到的，按公式算，信号的能量相等，按sum(x.^2)*T计算，10秒采集到的波形的能量更大。

现实情况中，比较两个波形的能量或有效值，都是采样率相同，采样时间相同，所有不会遇到如此纠结的问题。

生成一组信号：

fs=1000;

>> N=1000;

>> n=0:N-1;

>> t=n/fs;

>> x=sin(2*pi*100*t);

>> nfft=1024;

>> deltF=fs/nfft;

>> window=hanning(N);

>> %直接法，periodogram函数得到的功率谱密度

>>[Pxx_period,f_period]=periodogram(x,window,nfft,fs); 

> noverlap=50; 

>>[Pxx_welch,f_welch]=pwelch(x,window,noverlap,nfft,fs);

计算原始信号的有效值为：0.0224

画出频谱与功率谱密度为：



幅值谱的幅值理论上应为1，不到1的原因是fft变换的点数与采样点数不同所致。

利用FFT幅值谱的平方/N ,画功率谱密度结果跟上右图差不多。

xw=1.633*x.*window';     % 加汉宁窗(恢复系数为1.633)，能量修正系数使加窗后能量保证不变

mag=abs(fft(xw,nfft));

Pxx_1=mag.^2/N/fs;

f=(0:nfft/2-1)/nfft*fs;

plot(f,Pxx_11(1:512)*2),title('Pxx_11')

关于功率谱密度计算，先做自相关计算，再做FFT也能得到功率谱密度。

最后结果为：



总结：

当采样点数=nfft时，deltF*N/fs=1;

功率谱密度直接求和即是频域能量。

用幅值谱的平方估计频域能量时，除完点数，还要除以采样频率。

时域能量要*采样间隔（1/fs）

有效值的平方*采样时间=时域能量；

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频谱');




结果例：

形象一点说：亮度或灰度变化激烈的地方对应高频成分，如边缘；变化不大的地方对于低频成分，如大片色块区画个直方图，大块区域是低频，小块或离散的是高频把图像看成二维函数，变化剧烈的地方就对应高频，反之低频。

举个通俗易懂的例子：

一幅图象，你戴上眼镜，盯紧了一个地方看到的是高频分量

摘掉眼镜，眯起眼睛，模模糊糊看到的就是低频分量。  

图像的高低频是对图像各个位置之间强度变化的一种度量方法.

低频分量:主要对整副图像的强度的综合度量.

高频分量:主要是对图像边缘和轮廓的度量.

如果一副图像的各个位置的强度大小相等,则图像只存在低频分量,从图像的频谱图上看,只有一个主峰,且位于频率为零的位置.

如果一副图像的各个位置的强度变化剧烈,则图像不仅存在低频分量,同时也存在多种高频分量,从图像的频谱上看,不仅有一个主峰,同时也存在多个旁峰.

以上的现象可以通过对傅里叶变换的公式分析得出.

以下所说的积分是对x进行的.

exp(-jwx)的数值变化是均匀的,如果对exp(-jwx)进行积分,则积分值为零.如果对exp(-jwx)乘以一个加权函数f(x),则在对f(x)exp(-jwx)进行积分,积分值不一定为零.如果exp(-jwx)的取值为1时,则对f(x)exp(-jwx)积分,既为对f(x)积分,此时f(x)exp(-jwx)最大,既频谱中的主峰.如果f(x) 是常数则, 除w=0处f(x)exp(-jwx)的积分不为零外,在w不为零的其它处,f(x)exp(-jwx)的积分都为零. 

1.本文部分素材来源网络，版权归原作者所有，如涉及作品版权问题，请与我联系删除。

2.未经原作者允许不得转载本文内容，否则将视为侵权；

3.转载或者引用本文内容请注明来源及原作者；

4.对于不遵守此声明或者其他违法使用本文内容者，本人依法保留追究权等。

下面是我的个人微信公众号，关注【一个早起的程序员】精彩系列文章每天不断。



 

    
刚入手傅里叶变化，头大的很，尤其是经过傅里叶变化之后，看频谱，心里慌慌的，不过接触多了，也就明了的差不多了，回头在搜索资料，恨没有早点看到这个资料。
先说收获：
傅里叶变换之后，频谱有几个特点：
中心点是原图整幅图像的平均灰度
如果原图中有明显的横纹（竖纹），那么频谱图中就会有鲜明的竖线（横线）
看看几张图片：
通过控制傅里叶频谱中某些点，再观察变换回原图的状态，就能有一个比较好的理解了。



本文图来自：http://www.doc88.com/p-703756316355.html

转载于:https://www.cnblogs.com/OleNet/archive/2013/04/02/2996722.html