怎么加个窗函数啊
本帖最后由 kuai5237 于 2013-10-22 22:42 编辑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)) 怎么加个窗函数得到第二个图啊，大神们，帮忙改改程序啊
大神们拯救一下小弟啊 你说的是时域的窗函数，还是频率的窗函数。你确定你的画的功率谱而不是FFT谱？ sunyuxinhe 发表于 2013-10-25 14:02 static/image/common/back.gif
你说的是时域的窗函数，还是频率的窗函数。你确定你的画的功率谱而不是FFT谱？
是对测试的得到的加速度时程，经过fft变换，应该是功率谱吧，我也不是很明白啊，刚接触没多久。 加的频率的窗函数，你看第一图太毛躁了，需要加窗进行平滑处理。 窗函数就是对原信号 加权截断。你可以试试的。 与原始信号相乘 %% 窗函数测试
function main
clc
close all
Ts = 0.001;
Fs = 1/Ts;
%% 原始信号
t = 0:Ts:pi/2;
yt = sin(2*pi*5*t) + sin(2*pi*10*t) + sin(2*pi*15*t);
= Spectrum_Calc(yt, Fs);
figure
subplot(211)
plot(t, yt)
xlabel('t')
ylabel('y')
title('原始信号')
subplot(212)
plot(f, Yf)
xlabel('f')
ylabel('|Yf|')
xlim()
ylim()
title('原始信号频谱')
%% 加窗信号
win = hann(length(t));
yt1 = yt.*win';
= Spectrum_Calc(yt1, Fs);
figure
subplot(211)
plot(t, yt1)
xlabel('t')
ylabel('y')
title('加窗信号')
subplot(212)
plot(f1, 2*Yf1) % 2表示能量系数
xlabel('f')
ylabel('|Yf|')
xlim()
ylim()
title('加窗信号频谱')
end
%% 求取频谱
function = Spectrum_Calc(yt, Fs)
L = length(yt);
NFFT = 2^nextpow2(L);
Yf = fft(yt,NFFT)/L;
Yf = 2*abs(Yf(1:NFFT/2+1));
f = Fs/2*linspace(0,1,NFFT/2+1);
end
本帖最后由 hcharlie 于 2017-1-24 21:48 编辑
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早就离开论坛了，别的人知道道理就可以了。
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
页:
[1]