0
$\begingroup$

There is a data set which has two different sinusoids buried in noise. Data record has following properties:

A data records with a total data length of 2048 points. The sampling frequency is 100 Hz. Each record contains two sinusoids that are adjacent. Each sinusoidal signal is between 300 and 500 points long in a 2048 point record. Each sinusoid will have a different phase shifts and position in the record but will not overlap.

Solution: I have so far done 1. fft of dataset , 2. chopping the dataset and then apply fft , 3. average periodogram with 25% data overlap and 4.bandpass filter and selected (5-6) and (38.5-43.5) and selected two range range1 = 175:440; range2 = 1700:2048; but sill couldn't identify two frequency. Can anyone assist me whether my approach is correct or not?

%% Load the data
data = load('thouhid.txt');
fs = 100; % Sampling frequency
N = length(data);
time = (0:N-1) / fs;
% Perform FFT
fft_result = fft(data);
fft_magnitude = abs(fft_result);
frequencies = (0:length(data)-1) * fs / length(data);


 FFFT of the signal

window_size = 512;
overlap = 0;
step = window_size - overlap;

% Segment the data
windows = {};
for i = 1:step:N-window_size+1
    window = data(i:i+window_size-1);
    windows{end+1} = window;
end
f={};

ff={};
% Analyze each window
for i = 1:length(windows)
    window = windows{i};
    
    % Perform FFT
    fft_result = fft(window);
    fft_magnitude = abs(fft_result);
    frequencies = (0:length(window)-1) * fs / length(window);
    f{end+1}= frequencies;
    ff{end+1}=fft_magnitude;
    

end

 fft on chopped data

nfft = 256; % Number of FFT points, also the segment length
window_length = 256; % Length of each segment


% 25% overlap
noverlap_25 = round(0.25 * window_length);



% Calculate and plot the periodograms
figure;

pwelch(data, window_length, noverlap_25, nfft, fs);
title('Averaged Periodogram with 25% Overlap');

Average Periodogram

order = 4; % Filter order
filtered_signal_566 = bandpass_filter(data, 5, 6, fs, order);
filtered_signal_39 = bandpass_filter(data, 38.5, 43.5, fs, order);


% Define the specific ranges for each signal
range1 = 175:440;
range2 = 1700:2048;

function filtered_signal = bandpass_filter(data, lowcut, highcut, fs, order)
    
    [b, a] = butter(order, [lowcut highcut]/(fs/2), 'bandpass');
    
    filtered_signal = filtfilt(b, a, data);
end

Band Pass filtered signal

dataset: https://drive.google.com/file/d/1TNdrL-pXFKzSsC3wRSiSJN4R37mYXARF/view?usp=drive_link

$\endgroup$
8
  • $\begingroup$ I see that same question was asked by same user here dsp.stackexchange.com/questions/93470/… with lesser details. $\endgroup$
    – Arpit Jain
    Commented Apr 9 at 7:28
  • $\begingroup$ yes , deleted that question now, thanks $\endgroup$ Commented Apr 9 at 7:56
  • 2
    $\begingroup$ Can you post the data set (or a few sets). This is much easier to answer if you can play around with the actual data. It would also help if you state your requirements what precisely do you need to know frequency, amplitude, phase, and time window. There is a trade off between accuracy and effort $\endgroup$
    – Hilmar
    Commented Apr 9 at 12:26
  • $\begingroup$ check the dataset . I need to find out the frequency, amplitude, and phase of the two sinusoid thats in the dataset $\endgroup$ Commented Apr 9 at 19:34
  • $\begingroup$ Do you need to do this in an online or offline setting? Ie, is speed or accuracy more important? $\endgroup$
    – Baddioes
    Commented Apr 9 at 21:25

1 Answer 1

0
$\begingroup$

The problem that you are going to run into with this dataset is the extremely low SNR. Once you have obtained frequency estimates, the best way that you can find amplitude and phase estimates is by using the non-linear least squares (NLS) method under the assumption that you have the frequency estimates.

The problem then becomes what is the best way to find frequency estimates. You can use the NLS method, which can be solved by a maximum likelihood solver, such as expectation maximization, if the noise is white Gaussian noise. However, this is pretty involved, especially if the noise is not guaranteed to be white Gaussian. The unmodified periodogram provides approximations of the true frequency values with variances $\mathcal{O}(\frac{1}{N})$, if you assume that the frequency spacing between the two sinusoids is greater than $\frac{2\pi}{N}$ 1. Applying this more generally, we can assume that the peaks of any spectral estimate provide arbitrarily close estimates of the true frequency values of interest.

Interestingly, in this case, many of the different spectral estimates seem to conflict with one another since the SNR is so low. Most agree on a digital frequency of 0.391469, but the second one varies quite a bit. Here's a few spectral estimates I used Frequency Estimates via Spectral Estimate Peaks

My personal preference is MVM, but likely you will need to pick some combination of these techniques to determine the peaks to within a certain degree of confidence that you can live with. From there, you can use NLS to compute amplitude and phase estimates. This code is provided below by picking peaks from the periodogram.

% Approximate Nonlinear Least Squares Estimate of amplitude and phase
data = load('thouhid.txt');
fs = 100;
N = length(data);
time = (0:N-1) / fs;

Data = fftshift(fft(data));
frequencies = (-N/2:N/2-1)/N;

figure; plot(frequencies*fs,10*log10(Data.*conj(Data)/length(Data))); title('Periodogram')
[~,idx] = sort(20*log10(abs(reshape(Data,[],1))),1,'descend');
omegas_p = frequencies(idx(1:4));

B = exp(1i*(0:N-1).'*(omegas));

beta = B\data(:);
mag_est = abs(beta);
phase_est = angle(beta);

See 1 for more details on all of the methods used.

1 Spectral Analysis of Signals by Petre Stoica and Randolph Moses

$\endgroup$

Your Answer

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

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