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

4 次查看（过去 30 天）
Katarzyna Wieciorek
评论： Walter Roberson 2015-6-12
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?

请先登录，再回答此问题。

采纳的回答

Walter Roberson
Walter Roberson 2015-6-8
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.
  7 个评论
显示 5更早的评论
Walter Roberson
Walter Roberson 2015-6-12
One point in time is all you asked for. It finds the time at which the amplitude is greatest for that frequency.
Walter Roberson
Walter Roberson 2015-6-12
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.

请先登录，再进行评论。

更多回答（1 个）

Image Analyst
Image Analyst 2015-6-8
I think you need to look at the spectrogram() function in the Signal Processing Toolbox.
  3 个评论
显示 1更早的评论
Katarzyna Wieciorek
Katarzyna Wieciorek 2015-6-9
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');

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Time-Frequency Analysis 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by