这是我写的一个关于频谱细化的程序！希望和大家探讨！

hnyanhua

分辨率提高是因为频域范围缩小A倍，那么用实验来实现是话，是通过对滤波后的信号重采样来完成的吗？是用resample（f_2,1，fs/B）吗？这样分辨率才可以提高fs/B倍。
好像提示有错。
??? Error using ==> resample at 55
Q must be a positive integer.
.

songzy41

原帖由 hnyanhua 于 2008-6-11 10:00 发表
分辨率提高是因为频域范围缩小A倍，那么用实验来实现是话，是通过对滤波后的信号重采样来完成的吗？是用resample（f_2,1，fs/B）吗？这样分辨率才可以提高fs/B倍。
好像提示有错。
??? Error using ==> resample at 55
Q must be a positive integer.

分辨率提高，实际上是把采样频率减小了A倍，在同样的数据样点下，分辨率Δf=fs/N，fs减小了，N不变，Δf自然提高了。楼主在用resample（x,p,q）时，p和q都要为整数，所以fs/B可能不为整数造成了错误。

hnyanhua

原帖由 songzy41 于 2008-6-11 11:12 发表

分辨率提高，实际上是把采样频率减小了A倍，在同样的数据样点下，分辨率Δf=fs/N，fs减小了，N不变，Δf自然提高了。楼主在用resample（x,p,q）时，p和q都要为整数，所以fs/B可能不为整数造成了错误。

那么如果分辨率想提高A倍，就要使采样频率减小A倍，即此时信号采样频率为B，这样不会有混迭。
只是，这时需要由A来求出B，是这样的吗？由分辨率提高倍数推出滤波器带宽？
%我对23楼中的程序改进，具体如下：有什么不妥的地方请指出。呵呵谢谢了。
clear all;
close all;
N=input('取时间分隔的点数Ｎ＝');%N的取值需满足采样定律
t=linspace(0,1,N);dt=1/N-1;fs=N-1;%取信号长度Tp=1s，给出采样周期Ts=dt,采样频率fs=1/Ts=N-1
x=10*sin(2*pi*64*t)+10*sin(2*pi*250*t)+20*sin(2*pi*256*t)+30*sin(2*pi*258*t)+20*sin(2*pi*512*t);%信号最高频率fc=512HZ
figure(1)
stem(t,x);grid%画出原始信号
title('原始信号x(t)')
xlabel('t(s)');ylabel('x（t）');
f=linspace(0,N-1,N);
X=fft(x);
figure(2)
plot(f,abs(X)/max(abs(X))),grid
xlabel('f(HZ)');ylabel('|X(jf)|');
f0=input('需求的中心频谱f0=');
x1=x.*exp(-j*2*pi*t*f0);%'.*'为元素群运算
figure(3)
stem(t,x1);grid
title('频移后的信号x1(t)')
xlabel('t(s)');ylabel('x1(t)');
  
X1=fft(x1,N);
figure(4)
plot(f,abs(X1)/max(abs(X1))),grid
xlabel('f(HZ)');ylabel('|X1(jf1)|');
%axis([f0 fs 0 1]):'( 坐标选择有点问题，我想是显示【f0 fs】，但是显示图形不正确。不知道是什么原因。
B=input('理想底通滤波器的带宽B＝');
n=find(f<B/2);%求出底通滤波器带宽内的下标
f1=f(n);%取出中段频率
X2=X1(n);%取出中段频谱
figure(5)
plot(f1,abs(X2)/max(abs(X2))),grid%画出滤波后的频谱
xlabel('f1(HZ)');ylabel('|X2(jf1)|');
x2=ifft(X2,N);
figure(6)
stem(t,x2),grid%画出滤波后的波形
title('滤波后的信号')
xlabel('t(s)');ylabel('x2(t)');
y=resample(x2,1,20);% 对信号进行A倍的抽取
t1=resample(t,1,20);
figure(7)
stem(t1,y),grid%画出抽取后信号的波形
title('抽取后的信号y(t)')
xlabel('t1(s)');ylabel('y(t)');

Y=fft(y,N);% 对信号进行抽取后的 FFT 运算
figure(8)
plot(f,abs(Y)/max(abs(Y))),grid
xlabel('f(HZ)');ylabel('|Y(jf1)|');
%axis([f0 fs 0 1]):'( 这里也是坐标系选择有问题。取时间分隔的点数Ｎ＝1025需求的中心频谱f0=250理想底通滤波器的带宽B＝10

hnyanhua

figure(7)中
y=resample(x2,1,20);% 对信号进行A倍的抽取
t1=resample(t,1,20);
figure(7)
stem(t1,y),grid%画出抽取后信号的波形
title('抽取后的信号y(t)')
xlabel('t1(s)');ylabel('y(t)');
若对信号进行50倍抽取，原信号采样频率fs=1024，（因为fc=512）,带宽B=20.48,
这些参数的选择可以通过"输入时间分隔的点数Ｎ＝"来确定。是这样的？

songzy41

原帖由 hnyanhua 于 2008-6-11 14:30 发表

那么如果分辨率想提高A倍，就要使采样频率减小A倍，即此时信号采样频率为B，这样不会有混迭。
只是，这时需要由A来求出B，是这样的吗？由分辨率提高倍数推出滤波器带宽？
%我对23楼中的程序改进，具体如下：有什 ...

这两语句我都可以执行，不知楼主为什么执行不了？
axis([f0 fs 0 1]) %坐标选择有点问题，我想是显示【f0 fs】，但是显示图形不正确。不知道是什么原因。
axis([f0 fs 0 1]) %这里也是坐标系选择有问题。取时间分隔的点数Ｎ＝1025需求的中心频谱f0=250理想底通滤波器的带宽B＝10
得到的图如下，但是这两张图已没有什么意义了，在figure4中是已移频后的，f0 fs等于250 1024，不能画出所需要在0频率附近的图形；figure8中已下采样，但还用fs=1024，同时也是把所需要的信号移出在显示外。

songzy41

原帖由 hnyanhua 于 2008-6-11 14:56 发表
figure(7)中
y=resample(x2,1,20);% 对信号进行A倍的抽取
t1=resample(t,1,20);
figure(7)
stem(t1,y),grid%画出抽取后信号的波形
title('抽取后的信号y(t)')
xlabel('t1(s)');ylabel('y(t)');
若对信号进行5 ...

我想先不谈图7的问题，而说你的算法。x2是滤波后的数据，有1025点长，在20倍重采样后y只有52点长，然后补0，做1025点的FFT，得到图如下。可以看到有很多漪波，这些漪波的峰值是否算作ZFFT的信号呢？
在ZFFT中x2应长20500，下采样后使y长1025点，再做1025点的FFT，这样才是正确的ZFFT运算。

hnyanhua

原帖由 songzy41 于 2008-6-11 17:49 发表

这两语句我都可以执行，不知楼主为什么执行不了？
axis([f0 fs 0 1]) %坐标选择有点问题，我想是显示【f0 fs】，但是显示图形不正确。不知道是什么原因。
axis([f0 fs 0 1]) %这里也是坐标系选择有问题。取时间分 ...

是的，我明白了。
X的频谱显示【0,1024】,移频后，X1的频谱显示的是X频谱的【250，1024+250】，即X的250HZ移到了X1的零频，
我想的是，在X1的坐标系中显示出零频就是X频谱的f0=250HZ，不知道用什么方法实现？
同样figure8中也是这个意思，在figure8中的零频处显示X频谱的250HZ，使读者看到图形立即可以明白是对250HZ出的频谱进行了细化。
axis()语句不能实现，应该怎样做呢？不知道楼主明白了我的意思吗？

hnyanhua

原帖由 songzy41 于 2008-6-11 17:59 发表

我想先不谈图7的问题，而说你的算法。x2是滤波后的数据，有1025点长，在20倍重采样后y只有52点长，然后补0，做1025点的FFT，得到图如下。可以看到有很多漪波，这些漪波的峰值是否算作ZFFT的信号呢？
在ZFFT中x2应 ...

这些漪波的峰值当然不能算作ZFFT的信号，是我们不希望的。可以消除吗？
在ZFFT中x2应长20500，下采样后使y长1025点，再做1025点的FFT，这样才是正确的ZFFT运算。
按照楼主的意思，是该把x2=ifft(X2,N);改为x2=ifft(X2,20500);这样才是正确的ZFFT运算？
但是我在很多参考书上都看到对抽取后的信号补零，再作FFT运算，这又是怎么回事呢？

hnyanhua

我对上面的程序进行了改进，使分辨率提高fs/B倍，分辨率的选择是任意的吗？
显示图形有问题呀？
x2=ifft(X2,fs/B*N);t1=linspace(0,1,fs/B*N)
figure(6)
stem(t1,x2),grid%画出滤波后的波形
title('滤波后的信号')
xlabel('t1(s)');ylabel('x2(t1)');
y=resample(x2,1,fs/B);% 对信号进行fs/B倍的抽取
figure(7)
stem(t,y),grid%画出抽取后信号的波形
title('抽取后的信号y(t)')
xlabel('t(s)');ylabel('y(t)');

Y=fft(y,N);% 对信号进行抽取后的 FFT 运算
figure(8)
plot(f,abs(Y)/max(abs(Y))),grid
xlabel('f(HZ)');ylabel('|Y(jf)|');
%采样频率fs=1024需求的中心频谱f0=250理想底通滤波器的带宽B＝20.48

songzy41

原帖由 hnyanhua 于 2008-6-12 14:55 发表
我对上面的程序进行了改进，使分辨率提高fs/B倍，分辨率的选择是任意的吗？
显示图形有问题呀？
x2=ifft(X2,fs/B*N);t1=linspace(0,1,fs/B*N)
figure(6)
stem(t1,x2),grid%画出滤波后的波形
title('滤波后的信号 ...

楼主用了50倍的比例，但似乎有些错误。时间的设置：
t1=linspace(0,1,fs/B*N)
这样得dt1=B/(fs*N)，fs1=(fs*N)/B，即使经过B倍的下采样，比原来的fs还大。
我感觉楼主还没有完全了解ZFFT怎么去做，我在36层上提出x2应该是20481点（指20倍--36层上说是20500点，有误），这20481点不是靠补零得到的，应该是x本身就有20481点。下面我按这思路对程序作了较大的修改，及figure8的图。把该图与36层的图比较，可看到没有漪波，峰值明显，把2个个靠近的峰值分离了。
clear all; clc; close all;
N=input('取时间分隔的点数Ｎ＝');%N的取值需满足采样定律
t=linspace(0,1,N);dt=1/N-1;fs=N-1;%取信号长度Tp=1s，给出采样周期Ts=dt,采样频率fs=1/Ts=N-1
N1=20*(N-1)+1;
t1=(1:N1)/fs;
x=10*sin(2*pi*64*t1)+10*sin(2*pi*250*t1)+20*sin(2*pi*256*t1)+30*sin(2*pi*258*t1)+20*sin(2*pi*500*t1);%信号最高频率fc=512HZ
figure(1)
stem(t1,x);grid%画出原始信号
title('原始信号x(t)')
xlabel('t(s)');ylabel('x（t）');
f=linspace(0,N-1,N);
X=fft(x,N);
figure(2)
plot(f,abs(X)/max(abs(X))),grid
xlabel('f(HZ)');ylabel('|X(jf)|');
f0=input('需求的中心频谱f0=');
x1=x.*exp(-j*2*pi*t1*f0);%'.*'为元素群运算
figure(3)
stem(t1,x1);grid
title('频移后的信号x1(t)')
xlabel('t(s)');ylabel('x1(t)');
  
X1=fft(x1,N);
figure(4)
plot(f,abs(X1)/max(abs(X1))),grid
xlabel('f(HZ)');ylabel('|X1(jf1)|');
%axis([f0 fs 0 1]) %坐标选择有点问题，我想是显示【f0 fs】，但是显示图形不正确。不知道是什么原因。
XX1=fft(x1);
ff=(0:N1-1)*fs/N1;
B=input('理想低通滤波器的带宽B＝');
n=find(ff<B/2);%求出底通滤波器带宽内的下标
f1=f(n);%取出中段频率
X2=XX1(n);%取出中段频谱
x2=ifft(X2,N1);
y=resample(x2,1,20);% 对信号进行A倍的抽取
t2=resample(t1,1,20);
fs1=fs/20;
f3=(0:N-1)*fs1/N+250;
Y=fft(y,N);% 对信号进行抽取后的 FFT 运算
figure(8)
plot(f3,abs(Y)/max(abs(Y))),grid
xlabel('f(HZ)');ylabel('|Y(jf1)|');
set(gca, 'XTickMode', 'manual', 'XTick', [250, 256, 258, 260, 270, 280, 300]);
%输入参数分别为  1025  250  20.48

hnyanhua

原帖由 songzy41 于 2008-6-11 11:12 发表

分辨率提高，实际上是把采样频率减小了A倍，在同样的数据样点下，分辨率Δf=fs/N，fs减小了，N不变，Δf自然提高了。楼主在用resample（x,p,q）时，p和q都要为整数，所以fs/B可能不为整数造成了错误。

32楼中楼主提到分辨率提高是把采样频率减小了A（A=fs/B）倍，即分辨率Δf=fs/N（fs减小，N不变）提高A倍。
这样的话，可以通过设置fs，B，来确定分辨率提高倍数？
程序可以改写如下？
clear all; clc; close all;
N=input('取时间分隔的点数Ｎ＝');%N的取值需满足采样定律
f0=input('需求的中心频谱f0=');
B=input('理想低通滤波器的带宽B＝');
t=linspace(0,1,N);dt=1/N-1;fs=N-1;%取信号长度Tp=1s，给出采样周期Ts=dt,采样频率fs=1/Ts=N-1
%N1=20*(N-1)+1;改为
N1=（fs/B）*(N-1)+1;
t1=(1:N1)/fs;%我觉得是楼主这里增加了信号长度:@o ？
x=10*sin(2*pi*64*t1)+10*sin(2*pi*250*t1)+20*sin(2*pi*256*t1)+30*sin(2*pi*258*t1)+20*sin(2*pi*500*t1);%信号最高频率fc=512HZ
figure(1)
stem(t1,x);grid%画出原始信号
title('原始信号x(t1)')
xlabel('t1(s)');ylabel('x（t1）');
f=linspace(0,N-1,N);
X=fft(x,N);
figure(2)
plot(f,abs(X)/max(abs(X))),grid
xlabel('f(HZ)');ylabel('|X(jf)|');
%f0=input('需求的中心频谱f0=');
x1=x.*exp(-j*2*pi*t1*f0);%'.*'为元素群运算
figure(3)
stem(t1,x1);grid
title('频移后的信号x1(t1)')
xlabel('t1(s)');ylabel('x1(t1)');
  
X1=fft(x1,N);
f11=(0:N-1)*fs/N+250;
figure(4)
plot(f11,abs(X1)/max(abs(X1))),grid
xlabel('f11(HZ)');ylabel('|X1(jf11)|');
怎样用set（）语句设置坐标，使它显示效果和图8一样啊？

XX1=fft(x1);
ff=(0:N1-1)*fs/N1;%频域取样点数N1
%B=input('理想低通滤波器的带宽B＝');
n=find(ff<B/2);%求出底通滤波器带宽内的下标
f1=f(n);%取出中段频率
X2=XX1(n);%取出中段频谱
x2=ifft(X2,N1);
%y=resample(x2,1,20);% 对信号进行A倍的抽取
%t2=resample(t1,1,20);
%fs1=fs/20;改为
y=resample(x2,1,fs/B);% 对信号进行A倍的抽取
t2=resample(t1,1,fs/B);
fs1=B;
f3=(0:N-1)*fs1/N+250;
Y=fft(y,N);% 对信号进行抽取后的 FFT 运算
figure(8)
plot(f3,abs(Y)/max(abs(Y))),grid
xlabel('f(HZ)');ylabel('|Y(jf)|');
set(gca, 'XTickMode', 'manual', 'XTick', [250, 256, 258, 260, 270, 280, 300]);:@) 这条语句非常好，使显示出的图形非常清楚。
%输入参数分别为  1025  250  20.48

hnyanhua

原帖由 songzy41 于 2008-6-12 17:14 发表

楼主用了50倍的比例，但似乎有些错误。时间的设置：
t1=linspace(0,1,fs/B*N)
这样得dt1=B/(fs*N)，fs1=(fs*N)/B，即使经过B倍的下采样，比原来的fs还大。
我感觉楼主还没有完全了解ZFFT怎么去做，我在36层上提 ...

是的，我没有完全了解ZFFT怎么去做？
x2=ifft(X2,fs/B*N);t1=linspace(0,1,fs/B*N)
这里没有注意到fs1=(fs*N)/B增大，只是考虑了x2，t1点数相同。
呵呵，谢谢指出啊。
我怎么觉得楼主真的是增加了信号长度，虽然信号频谱图很清楚，但是时域图形显示不是很好？
呵呵，有解决方法吗？

songzy41

原帖由 hnyanhua 于 2008-6-12 20:15 发表
32楼中楼主提到分辨率提高是把采样频率减小了A（A=fs/B）倍，即分辨率Δf=fs/N（fs减小，N不变）提高A倍。
这样的话，可以通过设置fs，B，来确定分辨率提高倍数？

可以这样来确定分辨率提高倍数。

原帖由 hnyanhua 于 2008-6-12 20:15 发表
t1=(1:N1)/fs;%我觉得是楼主这里增加了信号长度:@o ？

ZFFT的实质是把信号长度增加了，它实际上提高分辨率是以信号长度为代价的。

原帖由 hnyanhua 于 2008-6-12 20:15 发表
怎样用set（）语句设置坐标，使它显示效果和图8一样啊？

用
set(gca, 'XTickMode', 'manual', 'XTick', [250, 400, 500, 524, 600, 774]);
因为分辨率不高，是无法分辨250,256,258。

原帖由 hnyanhua 于 2008-6-12 20:28 发表
但是时域图形显示不是很好？
呵呵，有解决方法吗？

当取B=20.48，x数组长51201，要显示清楚很难，但可以显示一部分，例如只显示1024个点。

hnyanhua

ZFFT的实质是把信号长度增加了，它实际上提高分辨率是以信号长度为代价的。:@o
真的是这样的？
看来我对ZFFT的实质没有弄清楚，有待学习，呵呵。
可以提供一些参考文献吗？呵呵，谢谢了。:loveliness:

songzy41

楼主可看一下王济和胡晓编 “MATLAB在振动信号处理中的应用” （中国水利水电出版社）一书，其中有一节是关于ZFFT，带有ZFFT的MATLAB程序。