将使用听诊器测量的心跳信号通过硬件(主要是放大器和截止频率为 100hz 的低通滤波器)传给计算机的声卡。现在用截止频率 100hz 过滤信号。下面给出了找到峰值和每分钟节拍的代码。代码仅适用于某些情况。请帮我找出错误
clear all
%input the signal into matlab
[x,fs]=wavread('heartbeat.wav');
figure(1)
subplot(2,1,1)
x1=x(:,2);
plot(x1(500:10000),'r-');
title('unfiltered input x(n),cut off frequency 0.0270,passband 60hz,stopband 70hz');
ylabel('amplitude in volts');
xlabel('number of samples')
grid on
%to filter the signal above 50-60 hz
order=4;
h=fir1(4,0.0270,hamming(order+1));
y=filter(h,1,x1);
subplot(2,1,2)
plot(y(500:10000),'b-')
title('filtered output y(n),cut off frequency 0.0270,passband 50hz,stopband 60hz');
ylabel('amplitude in volts');
xlabel('number of samples')
grid on
%sound(y,5000)
th = max(y) * 0.9; %So here I'm considering anything less than 90% of the max as not a real peak... this bit really depends on your logic of finding peaks though which you haven't explained
Yth = zeros(length(y), 1);
Yth(y > th) = y(y > th);
Ydiff = diff(Yth);
Ydiff_logical = Ydiff < 0;
Ypeaks = diff(Ydiff_logical) == 1;
p=sum(Ypeaks)
N = length(y);
duration_seconds=N/fs;
duration_minutes=duration_seconds/60;
BPM=p/duration_minutes;
bpm=ceil(BPM)
figure(2)
%frequency response of the filter
freqz(h,1)
title('Frequency response');
xlabel('normalized frequency (X pi) radians per sample');
ylabel('Magnitude');
grid on;