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.
.
songzy41
发表于 2008-6-11 11:12
本帖最后由 wdhd 于 2016-6-3 11:05 编辑
原帖由 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
发表于 2008-6-11 14:30
本帖最后由 wdhd 于 2016-6-3 11:05 编辑
原帖由 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=');%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】,但是显示图形不正确。不知道是什么原因。
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():'( 这里也是坐标系选择有问题。取时间分隔的点数N=1025需求的中心频谱f0=250理想底通滤波器的带宽B=10
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)');
若对信号进行50倍抽取,原信号采样频率fs=1024,(因为fc=512),带宽B=20.48,
这些参数的选择可以通过"输入时间分隔的点数N="来确定。是这样的?
songzy41
发表于 2008-6-11 17:49
本帖最后由 wdhd 于 2016-6-3 11:05 编辑
原帖由 hnyanhua 于 2008-6-11 14:30 发表
那么如果分辨率想提高A倍,就要使采样频率减小A倍,即此时信号采样频率为B,这样不会有混迭。
只是,这时需要由A来求出B,是这样的吗?由分辨率提高倍数推出滤波器带宽?
%我对23楼中的程序改进,具体如下:有什 ...
这两语句我都可以执行,不知楼主为什么执行不了?
axis() %坐标选择有点问题,我想是显示【f0 fs】,但是显示图形不正确。不知道是什么原因。
axis() %这里也是坐标系选择有问题。取时间分隔的点数N=1025需求的中心频谱f0=250理想底通滤波器的带宽B=10
得到的图如下,但是这两张图已没有什么意义了,在figure4中是已移频后的,f0 fs等于250 1024,不能画出所需要在0频率附近的图形;figure8中已下采样,但还用fs=1024,同时也是把所需要的信号移出在显示外。
songzy41
发表于 2008-6-11 17:59
本帖最后由 wdhd 于 2016-6-3 11:05 编辑
原帖由 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
发表于 2008-6-12 13:58
本帖最后由 wdhd 于 2016-6-3 11:05 编辑
原帖由 songzy41 于 2008-6-11 17:49 发表
这两语句我都可以执行,不知楼主为什么执行不了?
axis() %坐标选择有点问题,我想是显示【f0 fs】,但是显示图形不正确。不知道是什么原因。
axis() %这里也是坐标系选择有问题。取时间分 ...
是的,我明白了。
X的频谱显示【0,1024】,移频后,X1的频谱显示的是X频谱的【250,1024+250】,即X的250HZ移到了X1的零频,
我想的是,在X1的坐标系中显示出零频就是X频谱的f0=250HZ,不知道用什么方法实现?
同样figure8中也是这个意思,在figure8中的零频处显示X频谱的250HZ,使读者看到图形立即可以明白是对250HZ出的频谱进行了细化。
axis()语句不能实现,应该怎样做呢?不知道楼主明白了我的意思吗?
hnyanhua
发表于 2008-6-12 14:09
本帖最后由 wdhd 于 2016-6-3 11:05 编辑
原帖由 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
发表于 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('滤波后的信号')
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
发表于 2008-6-12 17:14
本帖最后由 wdhd 于 2016-6-3 11:06 编辑
原帖由 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=');%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】,但是显示图形不正确。不知道是什么原因。
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', );
%输入参数分别为102525020.48
hnyanhua
发表于 2008-6-12 20:15
本帖最后由 wdhd 于 2016-6-3 11:06 编辑
原帖由 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=');%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', );:@) 这条语句非常好,使显示出的图形非常清楚。
%输入参数分别为102525020.48
hnyanhua
发表于 2008-6-12 20:28
本帖最后由 wdhd 于 2016-6-3 11:06 编辑
原帖由 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
发表于 2008-6-13 08:20
本帖最后由 wdhd 于 2016-6-3 11:07 编辑
原帖由 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,256,258。
原帖由 hnyanhua 于 2008-6-12 20:28 发表
但是时域图形显示不是很好?
呵呵,有解决方法吗?
当取B=20.48,x数组长51201,要显示清楚很难,但可以显示一部分,例如只显示1024个点。
hnyanhua
发表于 2008-6-13 11:33
ZFFT的实质是把信号长度增加了,它实际上提高分辨率是以信号长度为代价的。:@o
真的是这样的?
看来我对ZFFT的实质没有弄清楚,有待学习,呵呵。
可以提供一些参考文献吗?呵呵,谢谢了。:loveliness:
songzy41
发表于 2008-6-13 18:22
楼主可看一下王济和胡晓编 “MATLAB在振动信号处理中的应用” (中国水利水电出版社)一书,其中有一节是关于ZFFT,带有ZFFT的MATLAB程序。