怎么加个窗函数啊

kuai5237

th=a1(1,1:75776);
f_th=fft(th,75776);
ff_th(1:3000)=0;
ff_th(3001:36372)=f_th(3001:36372);
ff_th(36372:75776)=0;
fth=ifft(ff_th,75776);
y=fth;
Y=fft(y,75776);
Pyy=Y.*conj(Y)/75776;
f=250/75776*(0:37888);
plot(f,Pyy(1:37889));
title('power spectral density')
xlabel('frequency(Hz)')
z=find(Pyy==max(Pyy))
fnew=f(z(1)) 怎么加个窗函数得到第二个图啊，大神们，帮忙改改程序啊
原始数据.rar (209.14 KB, 下载次数: 4)

2013-10-22 22:41 上传

点击文件名下载附件
数据

kuai5237

大神们  拯救一下小弟啊

sunyuxinhe

你说的是时域的窗函数，还是频率的窗函数。你确定你的画的功率谱而不是FFT谱？

kuai5237

sunyuxinhe 发表于 2013-10-25 14:02

                               
登录/注册后可看大图

你说的是时域的窗函数，还是频率的窗函数。你确定你的画的功率谱而不是FFT谱？

是对测试的得到的加速度时程，经过fft变换，应该是功率谱吧，我也不是很明白啊，刚接触没多久。   加的频率的窗函数，你看第一图太毛躁了，需要加窗进行平滑处理。

maieb2011

窗函数就是对原信号 加权截断。你可以试试的。

旧言虐心

与原始信号相乘

william

1. %% 窗函数测试
   
2. function main
   
3. clc
   
4. close all
   
5. Ts = 0.001;
   
6. Fs = 1/Ts;
   
7. %% 原始信号
   
8. t = 0:Ts:pi/2;
   
9. yt = sin(2*pi*5*t) + sin(2*pi*10*t) + sin(2*pi*15*t);
   
10. [Yf, f] = Spectrum_Calc(yt, Fs);
    
11. figure
    
12. subplot(211)
    
13. plot(t, yt)
    
14. xlabel('t')
    
15. ylabel('y')
    
16. title('原始信号')
    
17. subplot(212)
    
18. plot(f, Yf)
    
19. xlabel('f')
    
20. ylabel('|Yf|')
    
21. xlim([0 100])
    
22. ylim([0 1])
    
23. title('原始信号频谱')
    
24. %% 加窗信号
    
25. win = hann(length(t));
    
26. yt1 = yt.*win';
    
27. [Yf1, f1] = Spectrum_Calc(yt1, Fs);
    
28. figure
    
29. subplot(211)
    
30. plot(t, yt1)
    
31. xlabel('t')
    
32. ylabel('y')
    
33. title('加窗信号')
    
34. subplot(212)
    
35. plot(f1, 2*Yf1) % 2表示能量系数
    
36. xlabel('f')
    
37. ylabel('|Yf|')
    
38. xlim([0 100])
    
39. ylim([0 1])
    
40. title('加窗信号频谱')
    
41. end
    
42. %% 求取频谱
    
43. function [Yf, f] = Spectrum_Calc(yt, Fs)
    
44. L = length(yt);
    
45. NFFT = 2^nextpow2(L);
    
46. Yf = fft(yt,NFFT)/L;
    
47. Yf = 2*abs(Yf(1:NFFT/2+1));
    
48. f = Fs/2*linspace(0,1,NFFT/2+1);
    
49. end
    

复制代码

hcharlie

f_th=fft(th,75776);
ff_th(1:3000)=0;
ff_th(3001:36372)=f_th(3001:36372);
ff_th(36372:75776)=0;

你这几行有点问题，你的想法可能是FFT后低频清零，第1点是零频，除外是2：3000共2999点，把尾巴的2999点72778：75776也清零，如果也想高频清零，从36372：37888 对应的37890：39506，连同中点37889，共计36372：39506 清零就对了，你一塌刮子从36372到最后一起清零就错了，总之要适应FFT以后共轭对称的特点清零，否则IFFT以后会出现一大堆的虚部，以后大概就错了！！！
遗憾的是LZ早就离开论坛了，别的人知道道理就可以了。

zhangzy

hcharlie 发表于 2017-1-24 21:32
f_th=fft(th,75776);
ff_th(1:3000)=0;
ff_th(3001:36372)=f_th(3001:36372);

改成这样吗？   ff_th(36372：39506)=0