sampled signal fft signal

hi,i have a task that i dont get it how to be done.

i have sampled ecg signal with fs=1024 hz (vector )

the task is to plot the pulses per minute (1/min) of the ecg signal

s1 is the sampled vector

b=fft(s1)

a=length(s1)

w=0:(2*pi/a):2*pi*(a-1)/a)

well the axix on w,how the fs get in the work ,it must be that i need changing 1 min= 60s=60000ms

please help

i.e the signal is harmonic signal with doubles frequencies of the original frequency signal

[yuval ohayon]

anyone help/?

how do i scalinig the function

Answer 1

[David Goodmanson]

Hi yuval,

You don't have to bother with omega since the fft operates with straight frequency. I arbitrarily picked the number of points for a data set, and you supplied fs.

N = 10000;
fs = 1024;
delt = 1/fs;              % time array spacing, by definition
t = (0:N-1)*delt;

For the frequency array, the following is always true for ffts:

delf = fs/N;              % frequency array spacing (Hz), fft rule
f = (0:N-1)*delf;

Now for the fft and rescaled frequencies. The frequency spectrum plot is done in the usual way by taking absolute value, plotting positive frequencies only and doubling the result except fot the dc term. You will see peaks at 120 and its harmonics.

% first, make up some data, 120 bpm = 2 bps
f0 = 2;
y = exp(-10*mod(t,1/f0)+(.4)*rand(1,N));
figure(1)
plot(t,y)
% ----------------

% frequency domain calc
z = fft(y)/N;
f = f*60;                 % rescale frequencies from per sec to per min      <==

% plot 
halfN = floor((N+1)/2); 
f_plot = f(1:halfN);
absz_plot = abs(z(1:halfN));
absz_plot(2:end) = 2*absz_plot(2:end);      % don't double f = 0 (dc)
figure(2)
plot(f_plot,absz_plot)
xlim([0,1200])                              % optional

(Technically if N is even you should not double the last element in absz_plot either, but I will not hassle with that detail here. The whole process of half the frequencies, abs, double, etc. kind of sucks anyway).

[yuval ohayon]

david thank you,i havnt checked it yet but the answer seems very usefull ,god bless you, thanks

[yuval ohayon]

a=length(s1)

fs=1024

Tsample=1/1024 %%Ts=1/fs

t=(0:a-1)*Tsample %%time array

subplot(2,1,1)

plot(t,s1)

delf=fs/a %%%frequency array rule

f=(0:a-1)*delf ; %%'f array in sec'

F_w=fft(s1)/a

fs_sample_per_min=f*60 %rescale freq from per sec per min

b=floor(a+1)/2%%%%

b_axix=fs_sample_per_min(1:b)

Y_axix=abs(F_w(1:b))

subplot(2,1,2)

plot(b_axix,Y_axix,'b')

doesnt seem right ,the second graph should be with width ,david what wrong ?

[David Goodmanson]

I'm not sure what you mean by 'should be with width'. Hard to comment without seeing the s1 array. Could you attach that in a comment?

It's possible, though, that what you mean is you don't see much going on across the width of the frequency plot. For a typical ecg with a lot of pulses, most of the action is at the beginning of the frequency array. That's why I used the xlim command. Or if you have not done so, you could zoom in on the left region of your plot where frequencies of 0-1000 bpm are and see what is going on.

The b = ... command should be b = floor((a+1)/2) although what you have should only make a minor difference.

[yuval ohayon]

a=length(s1)

fs=1000

tsample=1/fs

t=(0:a-1)*tsample

figure(1)

plot(t,s1)

xlim([0 5])

xlabel('t[sec]')

ylabel('ecg(t) amplitude')

[yuval ohayon]

f=(0:a-1)*(fs/N) %%fft rule fs/n array

F_W=fft(s1)/N %%normelaize amplitude

F_per_min=f*60 %%rescale bps to bpm

b=floor(a+1)/2

b_axix=F_per_min(1:b)

Y_axix=abs(F_W(1:b))

figure(2)

plot(b_axix,Y_axix)

xlim([0 2000])

xlabel('bpm[rad/min]')

ylabel('F_W(ecg(w)')

[yuval ohayon]

david,why its neccesery changing a=length(s1)====>b,why i shorted the w axix by 2?

thats the diffference?

(if only need to change the f =60*f and normelaiz amplitude)