1
$\begingroup$

Recently, I came across a paper that focuses on signal processing using lock-in amplification.

the algorithm involves applying the Fast Fourier Transform (FFT) to a time-domain biosignal.

A band-pass filter is then applied in the frequency domain, specifically targeting the range of 0.7-1.5 Hz. This calculation is limited to spectral components with positive frequencies only.

The subsequent step involves performing the Inverse Fast Fourier Transform (IFFT) on the processed signal, resulting in a complex-numbered time-domain signal.

The article explains that the imaginary part of the signal is phase-shifted by π/2 relative to the real part of the signal.

Next step, the aurthor definiton a correlaion matrix is sum of the image pixel multipy the band-pass signal.Finally, use the dominate frequency reconstruct the signal like enter image description here

This can be used as a reference function for lock-in amplification.

While this approach is not very common, I'm having trouble understanding its implications. How can this be interpreted?

The band pass matlab code is

samplerate=30;      
for i = 1:1024
   filename = ['image' num2str(i) '.tiff'];
   img(:,:,i)= imread(filename);  %1024 points
   rawdata(i) = mean2(img(110:140,110:140,i))  % image ROI pixel mean as time signal
end
nfft=1024/2;
[row,column]=size(img(:,:,1));
ldf=samplerate/(nfft-1);          % calculate band pass location
Freq_range=[0.7 1.5];
R1 = round(Freq_range(1)/ldf);
R2 = round(Freq_range(2)/ldf);
Y=fft(rawdata);
Y(1:R1)=0;                       % band pass 
Y(R2:nfft+1)=0;
[m,p]=max(Y(1:nfft+1));          % the dominate frequency
Rcfeq_MXfeq=p*ldf;
BP_data=ifft(Y);                 % ifft signal is a complex number

%% lock in amplification
corrematrix = zeros(row,column);
for i = 1: length(rawdata)
    corrematrix = corrematrix +raw_img(:,:,i)*BP_data(i);
end

t = [0:1023]/samplerate; % 時間向量
for j = 1 : length(rawdata)
        H = real(corrematrix)*cos(2*pi*Rcfeq_MXfeq*t(j))+imag(corrematrix)*sin(2*pi*Rcfeq_MXfeq*t(j));

The artical is Photoplethysmographic imaging of high spatial resolution

DOI is 10.1364/BOE.2.000996

The first reply pic enter image description here

re the second reply pic enter image description here

Improve this question
$\endgroup$

1 Answer 1

Reset to default
1
$\begingroup$

The only positive frequencies came from the construction of the analytic signal. If we have a real signal (with positive and negative frequencies in its spectrum) $s(t)$ we can construct its analytic signal (complex signal and only positive frequencies) $s_a(t)$ by:

$s_a(t)=s(t)+j\hat{s}(t)$

where $\hat{s}(t)$ is the Hilbert transform of s(t). If you are not familiar with the Hilbert transform, I recommend you take a look at it, but basically, you are multiplying by $j$ the negative frequencies and by $-j$ the positive ones. You can take Fourier Transform on $s_a(t)$ and see that the negative frequencies are canceled

At this point you have an analytic signal. The analytic signal approach tends to be useful because it is easy to extract the instantaneous envelope and phase of the signal using complex number operations.

Improve this answer
$\endgroup$
5
  • $\begingroup$ It's helpful to know the problem is about Hilbert transform, but I'm not familiar with the Hilbert transform. Is the frequency in the reply meaning the real signal's FFT? ' you are multiplying by j the negative frequencies and by −j the positive ones. the FFT of the real signal? ' $\endgroup$
    – subborrrr
    Sep 9 at 10:28
  • $\begingroup$ I've tried But there is still valu in the negetive frequency fs = 30; % sampling frequency n=30 t = 0:(1/fs):(1-1/fs); % time vector signal = 3*cos(2*pi*5*t)+3*cos(2*pi*1*t); signal = signal+hilbert(signal)*1i; fftsig = fft(signal); siftfft = fftshift(fftsig); fshift = (-n/2:n/2-1)*(fs/n); % zero-centered frequency range powershift = abs(siftfft)*2/n; % zero-centered power figure plot(fshift,powershift) the result I put in the article $\endgroup$
    – subborrrr
    Sep 9 at 10:33
  • $\begingroup$ Could you tell me does my code is right to dealing with the band pass region as the paper . Or how to correct the code. Thanks! $\endgroup$
    – subborrrr
    Sep 9 at 10:40
  • $\begingroup$ In matlab y = hilbert(x) returns the analytic signal y of a real sequence x (not the Hilbert transform of x). You don't need signal = signal+hilbert(signal)*1i just signal_a=hilbert(signal). You can check the documentation Hilbert Transform hilbert $\endgroup$
    – Andy R
    Sep 9 at 13:21
  • $\begingroup$ I've tried signal_a=hilbert(signal) . The negative frequency is eliminated.But how could I comprehend the above result with the paper. Is it about the calculate the band pass in positive frequency? I'm stll comfuse. Thanks! $\endgroup$
    – subborrrr
    Sep 9 at 13:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.