时域插零的影响

1. 频域插值的实现 -- 序列补零后的FFT DFT末尾补零效果 原来x(n)的DFT为X(k)，在末尾补零rN个后成为y(n)，则 Y(k) = X(k/r)，在中间相当于插值，在频域插值。 然后，末尾补零x(n)的Z变换和DTFT是不会改变的仍为X(z), X(exp(jw)) % padding zero in the end, check effect to Xk clc; clear; f1=2.67; f2=3.75; f3=6.75; fs=20; N = 16; n = 0:N-1; x = sin(2*pi*f1*n/fs) + sin(2*pi*f2*n/fs) + sin(2*pi*f3*n/fs); X = fft(x, N); X = abs(X); k = fs/N * (0:N-1); %fs/N是频域分辨率 stem(k, X(1:N), 'o', 'r'); hold on % padding with rN zeros r = 2; x(N+1:r*N) = 0; X = fft(x, r*N); X = abs(X); k = fs/(r*N) * (0:r*N-1); stem(k, X(1:r*N), '*'); hold on % padding with rN zeros r = 4; x(N+1:r*N) = 0; X = fft(x, r*N); X = abs(X); k = fs/(r*N) * (0:r*N-1); stem(k, X(1:r*N), '>', 'y'); 2 时域零插值 -- 序列插值之后变化 % Interpolation 0, DFT effect clc; clear; f1=2.67; f2=3.75; fs=20; N = 16; n = 0:N-1; x = sin(2*pi*f1*n/fs) + sin(2*pi*f2*n/fs); X = fft(x, N); X = abs(X); k = fs/N * (0:N-1); %fs/N是频域分辨率 subplot(2,1,1); stem(k, X(1:N), 'o', 'r'); % 内插0，两个数之间插入1个0 xi = zeros(1, 2*N); % 看到没？MATLAB中，从大到小容易，从小到大，仍然是先构造大的 for i = 1:2:2*N xi(i) = x((i+1)/2); end X = fft(xi, 2*N); X = abs(X); k = fs/(2*N) * (0:2*N-1); %fs/N是频域分辨率 subplot(2,1,2);stem(k, X(1:2*N), 'o', 'r'); title('内插1个0后频谱'); 对采样信号的频谱,为提高计算效率,通常采用FFT算法进行计算,设数据点数为 　　N = T/dt = T.fs 　　则计算得到的离散频率点为 　　Xs(fi) , fi = i.fs/N , i = 0,1,2,…,N/2 　　这就相当于透过栅栏观赏风景,只能看到频谱的一部分,而其它频率点看不见,因此很可能使一部分有用的频率成分被漏掉,此种现象被称为栅栏效应. 不管是时域采样还是频域采样，都有相应的栅栏效应。只是当时域采样满足采样定理时，栅栏效应不会有什么影响。而频域采样的栅栏效应则影响很大，“挡住”或丢失的频率成分有可能是重要的或具有特征的成分，使信号处理失去意义。 减小栅栏效应可用提高采样间隔也就是频率分辨力的方法来解决。间隔小，频率分辨力高，被“挡住”或丢失的频率成分就会越少。但会增加采样点数，使计算工作量增加。解决此项矛盾可以采用如下方法：在满足采样定理的前提下，采用频率细化技术(ZOOM），亦可用把时域序列变换成频谱序列的方法。 例如：505Hz正弦波信号的频谱来说明栅栏效应所造成的频谱计算误差。 设定采样频率fs=5120Hz，软件中默认的FFT计算点数为512，其离散频率点为 fi = i.fs/N = i.5120/512=10×i ， i= 0，1，2，…，N/2 位于505Hz 位置的真实谱峰被挡住看不见，看见的只是它们在相邻频率500Hz或510Hz处能量泄漏的值。 若设 fs=2560Hz，则频率间隔df=5Hz，重复上述分析步骤，这时在505位置有谱线，我们就能得到它们的精确值。从时域看，这个条件相当于对信号进行整周期采样，实际中常用此方法来提高周期信号的频谱分析精度。 频谱泄露：截断信号时域上相当于是乘以了rectangular window，于是造成了频谱泄漏的问。 泄漏的原因来自两方面第一输入频率不是fs/n的整数倍，因为dft只能输出在fs/n的频率点上的功率，所以当输入频率不在fs/n的整数倍时，在dft的输出上就没有与输入频率相对应得点(dft输出是离散的)，那么输入频率就会泄漏到所有的输出点上，具体的泄漏分布取决于所采用的窗的连续域复利叶变换，对于没有使用窗的，相当于使用了矩形窗，矩形窗在进行连续傅立叶变换在一般的信号与系统书上都有。而对于非矩形窗，窗本身就会产生一定的泄漏，是通过加大主瓣的宽度来降低旁瓣的幅度，通常主瓣的宽度变成了矩形窗的两倍，例如当我们输入一个fs/n的整数倍的输入频率时，经过非矩形窗，dft输出会在两个fs/n的频点上有功率。 见参考书：lyon的understanding DSP. 旁瓣效应： 补零对频谱的影响： 进行zero padding只是增加了数据的长度，而不是原信号的长度。就好比本来信号是一个周期的余弦信号，如果又给它补了9个周期长度的0，那么信号并不是10个周期的余弦信号，而是一个周期的余弦加一串0，补的0并没有带来新的信息。其实zero padding等价于频域的sinc函数内插，而这个sinc函数的形状（主瓣宽度）是由补0前的信号长度决定的，补0的作用只是细化了这个sinc函数，并没有改变其主瓣宽度。而频率分辨率的含义是两个频率不同的信号在频率上可分，也就要求它们不能落到一个sinc函数的主瓣上。所以，如果待分析的两个信号频率接近，而时域长度又较短，那么在频域上它们就落在一个sinc主瓣内了，补再多的0也是无济于事的。 泄露是由于非整周期采样引起的，因为FFT最精确的是将周期信号映射到一个正交函数空间上（傅立叶变换常用三角函数空间），对周期性信号，只要是整周期采样（采样周期是信号周期的整数倍数描述不对，应该说采样时间长是信号周期的整数倍），是没有谱泄露的，对于非周期信号，无法达到整周期采样，所以总会有泄露，选择合适的窗函数可以控制泄露的严重程度。 混叠是采样频率与信号最高频率的关系引起的，满足采样定理，即采样频率〉=2倍信号最高频率，即可避免混叠，实际信号都受噪声干扰，白噪声是宽带的，所以采样频率即便很高，都不可避免地存在混叠，只是混叠程度小些，满足应用。再者，采样频率不是越高越好，我另外的帖子谈过，采样点数有限的情况下，采样频率与频率分辨率是相互矛盾的。 栅栏效应，频谱泄露，旁瓣效应 栅栏效应： 　　对采样信号的频谱,为提高计算效率,通常采用FFT算法进行计算,设数据点数为 　　N = T/dt = T.fs 　　则计算得到的离散频率点为 　　Xs(fi) , fi = i.fs/N , i = 0,1,2,…,N/2 　　这就相当于透过栅栏观赏风景,只能看到频谱的一部分,而其它频率点看不见,因此很可能使一部分有用的频率成分被漏掉,此种现象被称为栅栏效应. 不管是时域采样还是频域采样，都有相应的栅栏效应。只是当时域采样满足采样定理时，栅栏效应不会有什么影响。而频域采样的栅栏效应则影响很大，“挡住”或丢失的频率成分有可能是重要的或具有特征的成分，使信号处理失去意义。 减小栅栏效应可用提高采样间隔也就是频率分辨力的方法来解决。间隔小，频率分辨力高，被“挡住”或丢失的频率成分就会越少。但会增加采样点数，使计算工作量增加。解决此项矛盾可以采用如下方法：在满足采样定理的前提下，采用频率细化技术(ZOOM），亦可用把时域序列变换成频谱序列的方法。 例如：505Hz正弦波信号的频谱分析来说明栅栏效应所造成的频谱计算误差。 设定采样频率fs=5120Hz，软件中默认的FFT计算点数为512，其离散频率点为 fi = i.fs/N = i.5120/512=10×i ， i= 0，1，2，…，N/2 位于505Hz 位置的真实谱峰被挡住看不见，看见的只是它们在相邻频率500Hz或510Hz处能量泄漏的值。 若设 fs=2560Hz，则频率间隔df=5Hz，重复上述分析步骤，这时在505位置有谱线，我们就能得到它们的精确值。从时域看，这个条件相当于对信号进行整周期采样，实际中常用此方法来提高周期信号的频谱分析精度。 频谱泄露：截断信号时域上相当于是乘以了rectangular window，于是造成了频谱泄漏的问题。 在帖子上看到的解释： 泄漏的原因来自两方面第一输入频率不是fs/n的整数倍，因为dft只能输出在fs/n的频率点上的功率，所以当输入频率不在fs/n的整数倍时，在dft的输出上就没有与输入频率相对应得点(dft输出是离散的)，那么输入频率就会泄漏到所有的输出点上，具体的泄漏分布取决于所采用的窗的连续域复利叶变换，对于没有使用窗的，相当于使用了矩形窗，矩形窗在进行连续傅立叶变换在一般的信号与系统书上都有。而对于非矩形窗，窗本身就会产生一定的泄漏，是通过加大主瓣的宽度来降低旁瓣的幅度，通常主瓣的宽度变成了矩形窗的两倍，例如当我们输入一个fs/n的整数倍的输入频率时，经过非矩形窗，dft输出会在两个fs/n的频点上有功率。 见参考书：lyon的understanding DSP. 旁瓣效应： 补零对频谱的影响： 进行zero padding只是增加了数据的长度，而不是原信号的长度。就好比本来信号是一个周期的余弦信号，如果又给它补了9个周期长度的0，那么信号并不是10个周期的余弦信号，而是一个周期的余弦加一串0，补的0并没有带来新的信息。其实zero padding等价于频域的sinc函数内插，而这个sinc函数的形状（主瓣宽度）是由补0前的信号长度决定的，补0的作用只是细化了这个sinc函数，并没有改变其主瓣宽度。而频率分辨率的含义是两个频率不同的信号在频率上可分，也就要求它们不能落到一个sinc函数的主瓣上。所以，如果待分析的两个信号频率接近，而时域长度又较短，那么在频域上它们就落在一个sinc主瓣内了，补再多的0也是无济于事的。 连续时间信号经采样、截断后的序列为Xn(n)，其频谱函数XN(ejw)，并不随序列末端补零而改变，信号的频率分辨率为Fs/N.序列末端补零只能提高信号频谱显示的分辨率。换句话说，如果连续时间信号在离散化或时域加窗截断过程中，由于频谱泄漏或混叠等原因已造成信号频谱中信息的失真，则无论怎么补零做DFT，都无法再恢复已损失的信息。 　　提高信号的频率分辨率只有提高信号的采样频率或增加序列的截断长度N(信号的持续时间加长)。 1）数据后面补零-------不能提高信号的频率分辨率 　　序列末端补零后，尽管信号的频谱不会变化，但对序列做补零后L点DFT后，计算出的频谱实际上是原信号频谱在[0,2*pi)区间上L个等间隔采样，从而增加了对真实频谱采样的点数，并改变了采样点的位置，这将会显示出原信号频谱的更多的细节。故而数据后面补零可以克服栅栏效应。 2）数据间隔补零-------不能提高信号的频率分辨率 3）数据插值 相当于提高了信号的采样率，可以提高信号的频率分辨率 在数字信号处理中，高分辨率谱和高密度谱是较为易混淆的两个概念。获得高分辨率谱的途径是增加信号采样的时间tp，而高密度谱则是通过在时域补零得到的。高分辨率谱的用途很显示，可以分辨出频率间隔更小的两个频率分量，那么高密度谱有什么作用呢？要想明白高密度谱的概念，就不得知道一个名词：栅栏效应。高密度谱就是为了减小栅栏效效的。实际信号是无限长的，其频谱是连续的，但是要用计算机对信号进行频谱分析，就必须把它截短使之成为有限长度为tp的信号，这样的截短相当于对信号加矩形窗。经过加窗截取，信号的周期变为tp,其频谱相应地由原来的连续谱变为离散谱，离散谱的谱线只在f=1/tp的整数倍的位置上才出现，于是谱线间的实际信号的谱线有可能被挡住而损失掉，这称之为栅栏效应。例如截取信号长度为tp=0.5s，则可得到的谱线为2Hz，4Hz，6Hz，8Hz，…，若信号中包含频率为7Hz的分量，则该分量将被栅栏挡住，无法显示出来 （2*pi/T） 频率细化是70年代发展起来的一种新技术，其主要目的是识别谱图上的细微结构。从通常的FFT分析方法中我们已经知道，在频谱图上的有效频率分布范围是从0HZ到奈魁斯特频率fN为止，而谱线间隔(fs/N)决定了频率分辨能力，N表示数据点数，这里fs表示采样频率，且fN=fs/2。因此，要获得较高的分辨率可从下面两个方面进行。第一方面：降低采样频率，谱线间隔减小，但这样会降低奈魁斯特频率fN，从而导致频率分析范围小；第二方面：提高FFT计算长度N值，但这样要求较大的内存和降低运算速度。  高密度谱 不等于 高分辨率谱 数字域中，有一个概念是栅栏效应，通常可以靠补零的方式来减少此效应，但补零不能提高频率分辨率。 很多人在此很迷惑，在末尾加零后，使一个周期内的点数增加，必然使样点间隔更近，谱线更密，事以前看不到的谱分量就可以看到了，能够看到更多的谱，不是提高分辨力了吗？ 其实加零后，并没有改变原有记录的数据，原有数据的频谱一开始就存在，我们只是有的看不见，加零后只是让我们看见原来本就存在的频率，也就是说，原始数据代表的该有的频率就有，没有的频率加再多的零（极限是成连续的），也没法看见。 举个例子：用到MATLAB工具， 有信号x(n)=cos(0.48* pi*n)+cos(0.52*pi*n) 当0<=n<=19时，和0<=n<=99时，画出各自的信号图像 由上面的图可以知道，20点样本值后补零后频率成分没有变化，只是将没有显示的现实出来了，而取100个样本值后，可以清晰的看到信号中所含的两个频率分布在0.5 两边。20点的时候可以看到信号主要集中在0.5 处，但分辨不出0.48 和0.52 这两个具体值，20点补零后也不能分辨。 高频补零与过采样 IFFT在中间补的目的是提供频带的保护间隔，在时域产生的效果就是过采样。频率分辨率=1/时域持续时间长度时域采样率=1/频率长度 我们频域补零的前提是，频率分辨率不变。无论在高频或低频补零，都会使时域采样率增加，但根据乃奎斯特采样定理可知，频谱以采样率为间隔在频率上周期搬移，会使高频部分混跌，所以补零应在高频部分。 频域分辨率不变的高频域补零扩展，对应时域过采样（分辨率）补值。 clear all; clc; t=[0 1 2 3 4 5 6 7]; f=[0 1 2 3 4 5 6 7]; S=[1 2 0 4 3 5 2 6]; S1=[1 2 0 4 0 0 0 0 0 0 0 0 3 5 2 6]; Lt=length(t); Lf=length(S); Lf1=length(S1); Ts=1/Lf; Ts1=1/Lf1; figure(1); stem([-Lf/2:Lf/2-1],fftshift(S),'b+--'); hold on; f1=[0:Lf1-1]; stem(-Lf1/2:Lf1/2-1,fftshift(S1),'ro:'); hold off; figure(2); A=Lf*ifft(S); stem(t*Ts,A,'b+--'); hold on; t1=[0:Lf1-1]; A1=Lf1*ifft(S1); stem( t1*Ts1,A1,'ro:' ); hold off;