模仿频谱细化matlab程序

miao7mijao

模仿频谱细化：
大家帮我看看！

1. N=5120;
   
2. f=5120;
   
3. n=0:(N-1);
   
4. t=n/f;
   
5. x=sin(2*pi*50*t);
   
6. X=fft(x);
   
7. f=f*(1:N)/N;
   
8. subplot(2,1,1);
   
9. plot(f,abs(X),'-');
   
10. axis([0,115,0,3000]);
    
11. set(gca, 'XTickMode', 'manual', 'XTick', [0, 20,30,40,45,50,55,60,70,80,90,100,110,115]);
    
12. xlabel('Frequency');
    
13. ylabel('|F(k)|');
    
14. grid on;
    
15. 
    
16. f0=45;%移动的频率
    
17. h=exp(-j*2*pi*f0*t);
    
18. y=x.*h;%数字变频，把w=100点移动到原点
    
19. y1=resample(y,1,10);%重新采样，采样频率为fs/N
    
20. k1=0:1:511;
    
21. h1=exp((j*2*pi*f0*k1)/512);
    
22. y1=y1.*h1;%反移频，使得前后频率一致
    
23. y1=abs(fft(y1));%快速傅里叶变换
    
24. subplot(2,1,2);
    
25. plot(k1,y1);
    
26. axis([45,55,0,1000]);
    
27. set(gca, 'XTickMode', 'manual', 'XTick', [45,46,47,48,49,50,51,52,53,54,55]);
    
28. title('细化10倍后的频率特性');
    
29. xlabel('频率');
    
30. ylabel('幅度值');
    
31. grid on;

复制代码

为什么我细化之后的两个纵坐标的值不相等呢？两次最大值对应的横坐标值也不是同一个值？

songzy41

程序中有二处不合适:
1，f=f*(1:N)/N;不对，应从0开始；
2，二次FFT中，笫一次是5120点，笫二次是512点，点数不同，当然变换后幅值不相同。修改的办法是计算正弦信号的幅值，即在FFT后，乘以2/N。
修改后的程序是：

1. N=5120;
   
2. f=5120;
   
3. n=0:(N-1);
   
4. t=n/f;
   
5. x=sin(2*pi*50*t);
   
6. X=fft(x);
   
7. f=f*(0:N-1)/N;     % 被修改
   
8. subplot(2,1,1);
   
9. plot(f,abs(X)*2/N,'-');     % 被修改
   
10. axis([0,115,0,1.2]);     % 被修改
    
11. set(gca, 'XTickMode', 'manual', 'XTick', [0, 20,30,40,45,50,55,60,70,80,90,100,110,115]);
    
12. xlabel('Frequency');
    
13. ylabel('|F(k)|');
    
14. grid on;
    
15. 
    
16. f0=45;%移动的频率
    
17. h=exp(-j*2*pi*f0*t);
    
18. y=x.*h;%数字变频，把w=100点移动到原点
    
19. y1=resample(y,1,10);%重新采样，采样频率为fs/N
    
20. k1=0:1:511;
    
21. h1=exp((j*2*pi*f0*k1)/512);
    
22. y1=y1.*h1;%反移频，使得前后频率一致
    
23. y1=abs(fft(y1));%快速傅里叶变换
    
24. subplot(2,1,2);
    
25. plot(k1,y1*2/512);     % 被修改
    
26. axis([45,55,0,1.2]);     % 被修改
    
27. set(gca, 'XTickMode', 'manual', 'XTick', [45,46,47,48,49,50,51,52,53,54,55]);
    
28. title('细化10倍后的频率特性');
    
29. xlabel('频率');
    
30. ylabel('幅度值');
    
31. grid on;

复制代码

miao7mijao

谢谢你，楼上的！
我想知道为什么幅值要乘以2/N；
为什么幅值有那么大的变化，我的是2500左右，而改动之后却是1；

songzy41

可参看一下：http://forum.vibunion.com/thread-23665-1-1.html

miao7mijao

大家好！
我想问个问题，就是关于频率分辨率=fs/N和1/fs（我认为这个是采样数据的间隔），我觉的这两个是一个概念，有没有谁能给我解释一下呢？

songzy41

这是两个完全不同的概念，楼主可看一下，这两个量的单位也不一样，怎么会是“一个概念”呢？
fs/N常称作为频率分辨率，它实际是作FFT时谱图中的两条相邻谱线之间的频率间隔，也有称作步长。单位是Hz、Khz等。频率分辨率实际有二重含意，在这里只是其中一种。
1/fs的单位的s、ms、us或分、时...年等。1/fs代表采样周期，是时间域上两个相邻离散数据之间的时间差。
因此fs/N用在频率域，只在DFT以后的谱图中使用；而1/fs用时间域，只要数据经采样，离散化后任何其它的应用中都可使用。例如有的数字滤波器中就用到。

miao7mijao

真的太感谢你了！
讲的太好了
佩服！

微尼

请问songzy41，这个频谱细化程序在使用resample函数重新采样之前是不是没有滤波？

zhangnan3509

Δf=fs/N=1/T;Δf是频率采样间隔，同时也是频率分辨率的重要指标，如果这个值越小，则频率分辨率越高。
1/fs往往用在求时间序列上，如（0：N-1）*1/fs等等，如果这个不好理解，可以把前面的公式求倒数，这就清楚多了。
公式上的变换就可以把问题说得很清楚，而且不容易出错，这也是工科理论联系实际的魅力吧

[ 本帖最后由 zhangnan3509 于 2007-6-29 19:17 编辑 ]

zhangnan3509

你选的信号幅值是1，怎么会是2500？至于为什么乘以2/N
请看http://forum.vibunion.com/thread-23665-1-1.html
这里面songzy41 用公式推导的相当清楚，建议咱们再学习一下

wy558558558

这不就是加大频率分辨率了吗？:@o

zhangnan3509

首先，离散付立叶变换的定义本身比连续付立叶变换少了一个dt(采样时间间隔）；
然后，对于单频率成分的信号来说，经过矩形窗截断后的频谱在其信号频率处将放大T（做谱时间长度）倍，同样，对于相隔较远的多频率成分信号来说，相应的频率成分的幅值均将因截断而被放大T倍。
综合考虑这两种原因的话，也就是说我们用FFT做出的谱实际上是放大了T/dt=N（做谱点数）倍，因此，必须将此结果除以N。
（这是大版主yangzj的解释更加明了）这只不过是把该除的，除去，怎么会加大分辨率呢？打开上面的链接，看看推导结论，你就不会有这样的疑问了

[ 本帖最后由 zhangnan3509 于 2007-6-29 19:47 编辑 ]

wy558558558

不好意思就看图了，没仔细看内容:loveliness: 这和小波多分辩率很象．

zhangnan3509

频率细化方法提出的背景
频率细化是70年代发展起来的一种新技术，其主要目的是识别谱图上的细微结构。从通常的FFT分析方法中我们已经知道，在频谱图上的有效频率分布范围是从0HZ到奈魁斯特频率fN为止，而谱线间隔(fs/N)决定了频率分辨能力，N表示数据点数，这里fs表示采样频率，且fN=fs/2。因此，要获得较高的分辨率可从下面两个方面进行。第一方面：降低采样频率，谱线间隔减小，但这样会降低奈魁斯特频率fN，从而导致频率分析范围小；第二方面：提高FFT计算长度N值，但这样要求较大的内存和降低运算速度[8]。
在内存和FFT计算长度N有限制的情况下，既要不降低频率分析范围fN，而又要增加频率分辨率是矛盾的，为此出现了基于不同原理的各种选频细化分析方法，例如，扫频窄带分析法、基于复调制的ZFFT法、直接选抽法、级联FFT法、相位补偿细化和最大频谱的局部表示法等。最为常用的是复调制ZOOMFFT，相位补偿细化和级联三种方法。然而在计算效率、精度和灵活性等方面都比较理想的方法还是基于复调制的Zoom-FFT，因此得到了较多的应用。

zhangnan3509

1.复调制Zoom-FFT
复调制Zoom-FFT.输入信号为x(n),假设其频谱为|X(f)|,我们需要频率f0附近的频谱进行细微观察，则首先应对x(n)进行复调制，得到移频后的信号y(n),经过复调制后的信号y(n)的频谱是原来的频谱左移，欲观察的谱线已移至零频附近。这样就可以较低的频率对y(n)进行重新采样，为防止频谱混迭，在采样前应用理想低通滤波器进行滤波。
具体阐述如下：
（1）频移。为了将感兴趣的频段的下限频率移至原来的零频率位置，以便有可能将感兴趣频段放大到整个频率显示范围上，需首先对信号进行频率调制。这里采用的是复数调制法，如果欲将某一频率fo移至原来的零频处，则以原信号x1与 exp(-j2pi*f0*k*∆t) 相调制得：实部为 x1cos((2*pi*f0*k)/(N*∆f))，虚部为-x1sin(（2*pi*f0*k)/(N*∆f))。若令L0=f0/∆f(∆f---原有的频率分辨率)，即为频率在原频谱图中所对应的谱线序号，则实部和虚部即可以写为： x1cos(2*pi/N*L0*k)及 -x1sin(2*pi/N*L0*k)，合并实部和虚部可以得到调制后的信号为 wn=exp(-j*2pi/N)，

(2)滤波。数字低通滤波器是高截止特性的低通滤波器，可将从f0开始的一个所要求显示的窄频带f0到f1以外的所有频率成分滤掉，f0-f1仅为原截止频率的1/2^n(n=1,2,3,…)，此处2^n即为细化倍数,称之为细化因子。
(3)二次采样。二次采样是为了提高频率分辨率，使采样频率降至Fs/2^n(Fs是第一次采样的采样频率)。由采样定理可知，在采样个数仍为N时，采样频率下降为1/2^n，相当于总时间窗增长2^n倍，则频率分辨率亦将提高2n倍。这时的分辨率∆f’与原分辨率∆f之比为 1/2^n，经二次采样后的信号，进行复数FFT，便得到了细化的频谱。
由于细化2n倍时二次采样频率下降为原来的1/2^n，采样的记录长度亦应增至原来的2^n倍。应该指出，记录长度的增加仅在一次采样时增加了采样点数，而在完成二次采样后，点数仍为N，以后的FFT处理时间并未增加，因此，在细化2n倍时，计算时间并不会增加至2^n倍。
当然，移频法也有其缺点，就是一次分析仅能使指定的一段频谱得到细化与分析，而其余则均滤去，如欲进行多频段或全部频率范围内的细化，则要一次一次地进行重采样，然后再作预处理和分析，很费时间。

楼主用的应该是这种方法,resample有滤波的功能（谢谢songzy41先生)

[ 本帖最后由 zhangnan3509 于 2007-6-30 08:34 编辑 ]