Finding time points/intervals when frequency of the signal is about exact value

I would like to find at what time frequency of my signal is about 2 Hz. I know how to do fft of my signal but it does not let me find corresponding point in time. Do you have any idea how to do it?

Answer 1

[Walter Roberson]

http://www.mathworks.com/help/signal/ref/spectrogram.html

and look at Output Arguments.

Search the normalized frequencies output, w, for the one closest to your target frequency. That gives you a row offset into s. Search that row of s for the maximum absolute value. Index the column number into the time output, t, to get the time of the midpoint of that window.

At least that's what I figure from the documentation.

[Katarzyna Wieciorek]

Does it seem correct to you? Is there some simple way to avoid displaying the same values?

% Detrending signal of choosen electrode
sig_ori = val(el,:);
sig_det = detrend(sig_ori);

% Finding when frequency is about 2Hz
T = 60;                 
N = length(val);
fs = N/T;
NFFT = 2^nextpow2(length(sig_det));
[s_sig_det, w, ts] = spectrogram(sig_det,128,120,NFFT,fs,'yaxis');

[wrow,wcol,wv] = find(and(w>1.9,w<2.1));
[M,I] = max(abs(s_sig_det(wrow,:)));
t2Hz = ts(I);
h = msgbox(['Frequency was equal to 2Hz at time equal to ', num2str(t2Hz), ' s'],'Time');

[Katarzyna Wieciorek]

I used unique: t2Hz = unique(ts(I));

so the only problem is if it is correct.

[Walter Roberson]

targetfreq = 2.0
wrow = interp1(w, 1:length(w), targetfreq, 'nearest');

or

targetfreq = 2.0;
[~, wrow] = min(abs(w - targetfreq));

[Katarzyna Wieciorek]

Now it finds only one point in time, is it coincidence or code in this form finds only one?

Moreover, it's wrong :( it gives incorrect results.

[Katarzyna Wieciorek]

Any idea?

[Walter Roberson]

One point in time is all you asked for. It finds the time at which the amplitude is greatest for that frequency.

[Walter Roberson]

Are you looking for all the time points at which the largest peak is within a frequency range? How about situations where there is a "big" peak near the target frequency but the largest amplitude happens to be elsewhere, such as might occur if the peak is part-way between two bins so the energy is distributed to both of them.

Answer 2

[Image Analyst]

I think you need to look at the spectrogram() function in the Signal Processing Toolbox.

[Katarzyna Wieciorek]

I need something automatic. How can I extract wanted data from this matrix of complex doubles obtained? How its dimensions correspond to dimensions of my signal?

% Detrending signal of choosen electrode
sig_ori = val(el,:);
sig_det = detrend(sig_ori);

% Finding when frequency is about 2Hz
T = 60;                 
N = length(val);
fs = N/T;
NFFT = 2^nextpow2(length(sig_det));
s_sig_det = spectrogram(sig_det,128,120,NFFT,fs,'yaxis');

[Image Analyst]

Since you haven't uploaded "el", we can't really help you with your specific situation.

[Katarzyna Wieciorek]

All code:

%%PSIB 2015, Trabalho final, o problema 2
% Authors: Joao Leote, fc47337; Katarzyna Więciorek, fc45967

% This application adds a pulse with known temporal position, duration, 
% center frequency and bandwidth to a real signal downloaded from 
% the Physionet. STFT is used to make a time-frequency analysis and justify 
% all the chosen parameters. The limits of time and frequency location of 
% the selected method are explored by changing the pulse lengthor bandwidth.

%%Signal
% File with signals for each electrode
load S001R01_edfm.mat;

% Electrode selection
imshow('64electrodes.jpg');

prompt = 'Enter number of electrode. You can see them in separeted window.';
dlg_title = 'Electrode';
num_lines = 1;
def = {'11'};
answer = inputdlg(prompt,dlg_title,num_lines,def);
el = str2double(answer);

% If loop that checks if el is correct (in range from 1 to 64)
while(1)
    if (el >= 1) && (el <= 64)
        break;
    else
        prompt = 'Entered number is incorrect, type again:';
        dlg_title = 'Electrode';
        num_lines = 1;
        def = {'1'};
        answer = inputdlg(prompt,dlg_title,num_lines,def);
        el = str2double(answer);
    end
end

% Detrending signal of choosen electrode
sig_ori = val(el,:);
sig_det = detrend(sig_ori);

% Finding when frequency is about 2Hz
T = 60;                 
N = length(val);
fs = N/T;
NFFT = 2^nextpow2(length(sig_det));
[s_sig_det, w, ts] = spectrogram(sig_det,128,120,NFFT,fs,'yaxis');

%Searching the normalized frequencies output, w, for close to 2Hz 
targetfreq = 2.0;
wrow = interp1(w, 1:length(w), targetfreq, 'nearest');
%Searching found rows of spectogram for the maximum absolute value.
[M,I] = max(abs(s_sig_det(wrow,:)));
%Index the column number into the time output, t, to get the time of the midpoint of that window.
t2Hz = unique(ts(I));
h = msgbox(['Frequency was about 2Hz at time equal to: ', num2str(t2Hz),'seconds'],'Frequency 2Hz');
%%Sinusoidal pulse
% Sine
% Frequency
prompt = 'Enter frequency value: ';
dlg_title = 'Frequency';
num_lines = 1;
def = {'15'};
answer = inputdlg(prompt,dlg_title,num_lines,def);
freq = str2double(answer);

% Sine generation
T = 60;                 %Why 60? Why 5 in t2?
N = length(val);
fs = N/T;
t = linspace(0,T,T*fs);
t2 = linspace(0,5,5*fs);
sine = sin(2*pi*freq*t2);

%Pulse generation
pulsine = zeros(1,length(val));
pulsine (1,3201:4000) = 200*sine; %Why?

NFFT = 2^nextpow2(length(sine)); %http://www.mathworks.com/help/matlab/ref/nextpow2.html
sinefft = fft(sine,NFFT)/length(sine);    %Fast fourier transform
f = fs/2*linspace(0,1,NFFT/2+1);

% Signal + sinusoidal pulse
ssin = sig_det + pulsine;

%Graphics
figure(1)
subplot(2,1,1)
plot(t2,sine);
title('Sine')
xlabel('Time')
ylabel('Sine')
subplot(2,1,2)
plot(f,2*abs(sinefft(1:NFFT/2+1)))
title('Single-sided amplitude spectrum of sine')
xlabel('Frequency [Hz]')
ylabel('Amplitudes')

figure(2)
subplot(3,1,1)
plot(t,pulsine)
title('Sinusoidal pulse')
xlabel('Time')
ylabel('Sinusoidal pulse')

subplot(3,1,2)
plot(t,sig_det)
title('Detrended signal')
xlabel('Time')
ylabel('Signal')

subplot(3,1,3)
plot(t,ssin)
title('Signal with sinusoidal pulse added')
xlabel('Time')
ylabel('Signal + pulse')

figure(3)
spectrogram(ssin,128,120,NFFT,fs,'yaxis');

%http://www.mathworks.com/help/signal/ref/spectrogram.html

%%Gaussian-modulated sinusoidal pulse
tc = gauspuls('cutoff',freq,0.1,[],-40); %http://www.mathworks.com/help/signal/ref/gauspuls.html
tct = -tc : 1/fs : tc; %What is this?
t1 = linspace(-2.5,2.5,5*fs); %Why?
pg = gauspuls(t1,freq,1/freq); 
pulsga = zeros(1,length(val));
pulsga(1,3201:(3201+length(pg)-1)) = 200*pg; %Why?

NFFT = 2^nextpow2(length(pg));
gaussfft = fft(pg,NFFT)/length(pg);
f = fs/2*linspace(0,1,NFFT/2+1);

% Signal + Gaussian-modulated sinusoidal pulse
spg = sig_det+pulsga;

figure(4)
subplot(2,1,1)
plot(t1,pg);
title('Gaussian-modulated sinusoidal pulse')
xlabel('Time')
ylabel('Amplitude')

subplot(2,1,2)
plot(f,2*abs(gaussfft(1:NFFT/2+1)));
title('Single-sided amplitude spectrum of Gaussian-modulated sinusoidal pulse')
xlabel('Frequency [Hz]')
ylabel('Amplitude')

figure(5)
subplot(3,1,1)
plot(t,pulsga)
xlabel('Time'); 
ylabel('Pulse');
title('Gaussian-modulated sinusoidal pulse');

subplot(3,1,2)
plot(t,sig_det)
xlabel('Time'); 
ylabel('Signal');
title('Detrended signal');

subplot(3,1,3)
plot(t,spg)
xlabel('Time'); 
ylabel('Signal + pulse');
title('Signal with Gaussian-modulated sinusoidal pulse added');

figure(6)
spectrogram(spg,128,120,NFFT,fs,'yaxis');