刚入手傅里叶变化,头大的很,尤其是经过傅里叶变化之后,看频谱,心里慌慌的,不过接触多了,也就明了的差不多了,回头在搜索资料,恨没有早点看到这个资料。
先说收获:
傅里叶变换之后,频谱有几个特点:
- 中心点是原图整幅图像的平均灰度
- 如果原图中有明显的横纹(竖纹),那么频谱图中就会有鲜明的竖线(横线)
看看几张图片:
通过控制傅里叶频谱中某些点,再观察变换回原图的状态,就能有一个比较好的理解了。
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.对于不遵守此声明或者其他违法使用本文内容者,本人依法保留追究权等。
下面是我的个人微信公众号,关注【一个早起的程序员】精彩系列文章每天不断。
刚入手傅里叶变化,头大的很,尤其是经过傅里叶变化之后,看频谱,心里慌慌的,不过接触多了,也就明了的差不多了,回头在搜索资料,恨没有早点看到这个资料。
先说收获:
傅里叶变换之后,频谱有几个特点:
- 中心点是原图整幅图像的平均灰度
- 如果原图中有明显的横纹(竖纹),那么频谱图中就会有鲜明的竖线(横线)
看看几张图片:
通过控制傅里叶频谱中某些点,再观察变换回原图的状态,就能有一个比较好的理解了。
转载于:https://www.cnblogs.com/OleNet/archive/2013/04/02/2996722.html