|
发表于 2023-10-13 19:51:49
|
显示全部楼层
为y函数加上汉宁窗,可以使用以下代码:
w = hann(N)';
y = y .* w;
这里使用了hann函数生成长度为N的汉宁窗,并将其与y相乘,实现了对y的汉宁窗校正。
数字滤波可以使用fir1函数实现,代码如下:
fc = 15;
[b, a] = fir1(100, fc/(fs/2));
y = filter(b, a, y);
这里使用fir1函数生成了一个长度为101的FIR低通滤波器,截止频率为15Hz,然后使用filter函数对y进行滤波。
时频分析可以使用短时傅里叶变换(STFT)实现,代码如下:
window = hann(256)';
noverlap = 128;
nfft = 512;
[S, f, t] = spectrogram(y, window, noverlap, nfft, fs);
S = abs(S);
S = 20*log10(S/max(S(:)));
这里使用了spectrogram函数对y进行STFT分析,生成了时频图S,其中f是频率轴,t是时间轴。最后将S进行归一化并转化为分贝单位。
完整代码如下:
fs=180;
N=1024;
n=0:N-1;
t=n/fs;
y=0.68*sin(2*pi*10*t+pi/2)+0.45*cos(2*pi*20*t+1/4*pi)+0.28*randn(size(t));
w = hann(N)';
y = y .* w;
fc = 15;
[b, a] = fir1(100, fc/(fs/2));
y = filter(b, a, y);
window = hann(256)';
noverlap = 128;
nfft = 512;
[S, f, t] = spectrogram(y, window, noverlap, nfft, fs);
S = abs(S);
S = 20*log10(S/max(S(:)));
figure(1)
subplot(3,1,1)
plot(t,y)
xlabel('time/s')
ylabel('amplitude/1')
axis([0,1,-1.5,2])
title('0.68*sin(2*pi*10*t+pi/2)+0.45*cos(2*pi*20*t+1/4*pi)+randn')
a=fft(y,N);
b=abs(a/N);
b1=b(1:N/2+1);
b1(2:end-1)=2*b1(2:end-1);
f=(0:(N/2))*fs/N;
subplot(3,1,2)
plot(f,b1)
xlabel('frequency/Hz')
ylabel('amplitude/1')
title('A 0.68 0.45 fs 10 20')
axis([0,50,0,1])
subplot(3,1,3)
imagesc(t, f, S)
axis xy
xlabel('time/s')
ylabel('frequency/Hz')
title('spectrogram')
colorbar |
|