声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

12
返回列表 发新帖
楼主: Jnny_CN

[FFT] [求助]细化一组含有谐波的正弦波数据

[复制链接]
 楼主| 发表于 2008-6-2 12:33 | 显示全部楼层
本帖最后由 VibInfo 于 2016-10-13 16:47 编辑
原帖由 zhwang554 于 2008-6-1 17:16 发表
如果用离散频谱校正法(最近论坛贴子上介绍能量中心法,比值法程序)求出基波及其4个谐波的振幅值,相位值,频率值,重构一个信号,这个信号和原信号基本一样,但振幅谱的谱波清晰,


下图是500.0txt结果
兰色原信号的振 ...

不知您可否把具体的过程发表出来,让我学习一下,看了您的结果真是不错。谢谢。
回复 支持 反对
分享到:

使用道具 举报

发表于 2008-6-3 05:31 | 显示全部楼层
请问在上面程序段中“fre1=(fre*k-300)/(fs/2);”其中fre代表的是什么,程序中怎么没有体现,我试着仿真找不到fre,请楼主说明,谢谢。
发表于 2008-6-3 08:39 | 显示全部楼层
fre在给出的程序中没有表明,它和输入的文件名有关,在上例子中fre=1000。
发表于 2008-6-5 20:04 | 显示全部楼层

那个细化频谱是怎么一回事?有具体的实现程序吗?

最后对经过滤波器的信号抽取,进行快速傅里叶变换,傅里叶变换的点数是怎么选择的呢?
假如有这样的程序:
clear all;
close all;
N=input('取时间分隔的点数N=');
t=linspace(0,1,N);dt=1/N;%信号长度1秒,给出时间分割
f=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);
figure(1)
subplot(2,1,1);
stem(t,f,'g');grid%画出原始信号
title('原始信号')
xlabel('t(s)');ylabel('ft');
fw=input('需求的频谱宽度fw=');
Nf=input('需求的频谱点数Nf=');
f1=linspace(0,fw,Nf);df=fw/(Nf-1);dw=2*pi*df;
F=f*exp(-j*t'*2*pi*f1)*dt;%f傅立叶变换
subplot(2,1,2);
plot(f1,abs(F)/max(abs(F))),grid
xlabel('f1(HZ)');ylabel('|F(jf1)|');


f_0=input('需求的中心频谱f_0=');
f_1=f.*exp(-j*2*pi*t*f_0);%'.*'为元素群运算
figure(2)
subplot(2,1,1);
stem(f_1);grid
title('频移后的信号')
xlabel('t(s)');ylabel('f_1(t)');
F1=f_1*exp(-j*t'*2*pi*f1)*dt;%f_1傅立叶变换
subplot(2,1,2);
plot(f1,abs(F1)/max(abs(F1))),grid
xlabel('f1(HZ)');ylabel('|F1(jf1)|');



fw0=input('理想底通滤波器的带宽fw0(HZ)');
n2=find((f1>-fw0)&(f1<fw0));%求出底通滤波器带宽内的下标
f2=f1(n2);%取出中段频率
F2=F1(n2);%取出中段频谱
figure(3)
subplot(2,1,2);
plot(f2,abs(F2)/max(abs(F2))),grid%画出滤波后的频谱
xlabel('f2(HZ)');ylabel('|F2(jf2)|');
f_2=F2*exp(j*2*pi*f2'*t)/pi*dw;%对中段频谱求傅立叶逆变换
subplot(2,1,1);
stem(f_2),grid%画出滤波后的波形
title('滤波后的信号')
xlabel('t(s)');ylabel('f_2(t)');



y=decimate(f_2,20);% 对信号先进行滤波,然后进行 20 倍的抽取
figure(4)
subplot(2,1,1);
stem(y),grid%画出滤波后的波形
title('抽取后的信号')
xlabel('t(s)');ylabel('y(t)');
Y=fft(y);% 对信号进行 20 倍抽取后的 FFT 运算
subplot(2,1,2);
plot(abs(Y)/max(abs(Y))),grid
xlabel('f1(HZ)');ylabel('|Y(jf1)|');









程序和图形.doc

125.5 KB, 下载次数: 22

发表于 2008-6-5 20:09 | 显示全部楼层

那个细化频谱是怎么一回事?有具体的实现程序吗?

是对进行 20 倍抽取后的信号补零进行快速变换吗?
发表于 2008-6-14 10:35 | 显示全部楼层

回复16楼

用校正法求出各谐波的振幅值,相位值和频率值, 重构原信号,其振幅谱干净了,噪音去除了, 但这种方法分析谐波成份仍不能消除泄漏的影响,重构时泄漏也包含在其中了.
    正确的方法还是设法降低泄漏的影响,如加长数据,合适的取样频率等.
   
    不过这是一种滤波方法, 一种不用滤波器的去噪方法. 用校正方法求出一个信号的各种频率成分的振幅值,相位值和频率值, 重构原信号, 噪音去除了, 而不必设计多个窄带滤波器滤出各频率成份组成原信号..

    您希望把具体的过程发表出来,学习一下。

   程序如下,可能有错,谨供参考

close all;clc;clear all;
N=2048;
yy=load('500.0.txt');
y=yy'-mean(yy);
y1 = y(N:2*N-1);
win =  hanning(N)';;
win1 = win/sum(win);
y11= y1.*win1;
y11_fft = fft(y11,N);
a1 = abs(y11_fft);
p1 = mod(phase(y11_fft),2*pi);
y2 = y(1:2*N-1);
win =  hanning(N)';
winn =  conv(win,win);%apFFT须要卷积窗
win2 = winn/sum(winn);
y22= y2.*win2;
y222=y22(N:end)+[0 y22(1:N-1)];%构成长N的apFFT输入数据
y2_fft = fft(y222,N);
a2 = abs(y2_fft);
p2=mod( phase(y2_fft),2*pi);%相位值,不必校正
     ee=(p1-p2)/pi/(1-1/N);%频率校正
      aa=(a1.^2)./a2;%振幅校正
plot(10*log10(a1(2:N/16)),'b');
hold
      index1=16; index2=32;   index3=48;   index4=63;   index5=79;%从振幅谱比找出各谐波峰值处的频率
     
    disp('相位校正值')
       p11=p2(index1+1)
       p22=p2(index2+1)
       p33=p2(index3+1)
       p44=p2(index4+1)
       p55=p2(index5+1)
    disp('频率校正值')
       f1=(index1+ee(index1+1))
       f2=(index2+ee(index2+1))
       f3=(index3+ee(index3+1))
       f4=(index4+ee(index4+1))
       f5=(index5+ee(index5+1))
    disp('振幅校正值')
       z1=aa(index1+1)
       z2=aa(index2+1)
       z3=aa(index3+1)
       z4=aa(index4+1)
       z5=aa(index5+1)
    t=(-N+1:2*N-1)/N;
    s2=(z1*cos(2*pi*f1*t+p11)+z2*cos(2*pi*f2*t+p22)+z3*cos(2*pi*f3*t+p33)+z4*cos(2*pi*f4*t+p44)+z5*cos(2*pi*f5*t+p55))*2;%重构信号
Sin3=s2(1:2*N-1);
Out3=apFFT(Sin3,N);
a3=abs(Out3);
plot(10*log10(a3(2:N/16)),'r');
title('500.0txt的对数振幅谱  FFT 兰  用FFT和apFFT组合校正法重组 红');
set(gca,'XTick',0:f1:N/16)
set(gca,'XTickLabel',{0:f1*64710/N:4500})
ylim([-150,0]);
grid

[ 本帖最后由 zhwang554 于 2008-6-14 15:02 编辑 ]
 楼主| 发表于 2008-6-16 16:40 | 显示全部楼层
原帖由 zhwang554 于 2008-6-14 10:35 发表
用校正法求出各谐波的振幅值,相位值和频率值, 重构原信号,其振幅谱干净了,噪音去除了, 但这种方法分析谐波成份仍不能消除泄漏的影响,重构时泄漏也包含在其中了.
    正确的方法还是设法降低泄漏的影响,如加长数据,合 ...


谢谢您的帮助,令我受益匪浅:handshake
发表于 2011-4-1 11:09 | 显示全部楼层
回复 21 # zhwang554 的帖子

请问这个数的点数是怎么取的?N为什么是2048?它的数据共7999个点,其中一个周期128个点,取2048是为了什么?为了全相位周期延拓么?
发表于 2011-4-1 21:58 | 显示全部楼层
本帖最后由 zhwang554 于 2011-4-1 21:58 编辑

回复 23 # qiuyun0214 的帖子


N阶全相位fft需2N-1个数据, 现数据共7999个点,只能作4000阶apfft, 若为了有快速DFT, 只能取2048阶,再上一阶FFT是4096,,需要8096个数据, 超过数据7999个点
全相位fft数据不能通过周期延拓得到,必连续采样得到
发表于 2011-4-6 13:46 | 显示全部楼层
回复 24 # zhwang554 的帖子

谢谢老师!明白这个数据是怎么取的了!
发表于 2011-5-15 11:14 | 显示全部楼层
xialaixuexi
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-18 10:10 , Processed in 0.066270 second(s), 20 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表