convert the time domain signal into frequency domain signal

Hi, I have represented the acceleration data of 3 axes(x, y, and z) in time domain as shown in the Graph

I would like to extract from the acceleration data some measurements (e.g. mean, energy, entropy and correlation) in the frequency domain. Therefore, I applied FFT transform in order to convert the time domain signal into frequency domain signal.

    xx=  fft(x_Segments{1});
    plot(xx)
    yy=  fft(y_Segments{1});
    plot(yy)
     zz=  fft(z_Segments{1});
     plot(zz)

However, the resulted graphs make no sense( so strange) which is not expected at all !! This is an example of the frequency domain signal of the x-axis.

Note, x_Segments, y_Segments, and z_Segments contain the data of X, Y, Z axes respectively

As we can see the graph doesn't make sense, so, is there any steps that I should follow before using the FFT function to get the expected frequency domain signal?

I really appreciate your help guys. thanks

Answer 1

[Star Strider]

See this version of the documentation for fft.

[neamah al-naffakh]

Hi mate, I know that and already use it (you can see that in my question) but the graph that I got is not what I'm expecting. Regards.

[Star Strider]

You’re not using it correctly.

The reason I directed you to that particular version of the fft documentation is that it will show you how to produce the correct plot. That will be what you’re expecting.

[neamah al-naffakh]

Dear Star Strider, I have used this one

       Y = fft(x)

which is equivalent to this

     xx=  fft(x_Segments{1});

x_Segments{1} is my time domain signal and xx is the frequency domain shall I use another one? Thanks for your help. Regards.

[Star Strider]

My pleasure.

*Please read the documentation for the fft* function that I linked to in my Answer (and again here). Specifically see the code between the first (top) and second plot figures. It will tell you all you need to know about plotting and interpreting the output of your fft call.

You have to plot the absolute value in order to make any sense of it, and select the appropriate frequency range (specifically 0 Hz up to the Nyquist frequency) in your plot.

[neamah al-naffakh]

I'm really appreciate your help bellow is the step that i did with the following graphs, please let me know if there is any comments or correction from your side

    plot( X_Segments{1})
title('Noisy signal');
xlabel('Samples');
ylabel('Amplitude');

     %plot magnitude spectrum of the signal
X_mags = abs(fft(X_Segments{1}));
figure(2)
plot(X_mags)
xlabel('DFT Bins');
ylabel('Magnitude')

    figure(3)
num_bins = length(X_mags);
plot([0:1/(num_bins/2 -1):1], X_mags(1:num_bins/2))
xlabel('Normalised frequency (\pi rads/sample)')
ylabel('Magnitude')

regards

[Star Strider]

Your figure(3) looks essentially correct to me. I have only the plot of your original time-domain signal, not the signal itself.

The only change I would make is to multiply ‘X_mags’ by 2/length(X_Segments{1}) to normalise for the one-sided amplitude, and energy in the signal.

[neamah al-naffakh]

Dear Star Strider , thanks, it's really so kind of you. do you mean like this? if no , could you please make the correction for me because I'm not quite syre what do you mean?

 plot( X_Segments{1})
title('Noisy signal');
xlabel('Samples');
ylabel('Amplitude');

then

   %plot magnitude spectrum of the signal
X_mags = abs(fft(X_Segments{1}));
figure(2)
plot(X_mags)
xlabel('DFT Bins');
ylabel('Magnitude')
 X_mags=X_mags*(2/length(X_Segments{1}));

 figure(3)
num_bins = length(X_mags);
plot([0:1/(num_bins/2 -1):1], X_mags(1:num_bins/2))
xlabel('Normalised frequency (\pi rads/sample)')
ylabel('Magnitude')

[Star Strider]

My pleasure.

The change I would do in figure(3) is:

plot([0:1/(num_bins/2 -1):1], X_mags(1:num_bins/2)*2/length(X_Segments{1}))

[neamah al-naffakh]

dear Star Strider, you are so helpful, thank you sir . Sir, now the next step is to measure (mean, energy, entropy and correlation) in the frequency domain, shall I use this matrix (X_mags), which I have calculated in this function

     X_mags = abs(fft(X_Segments{1}));

or i should use this one?

     X_mags=X_mags*(2/length(X_Segments{1}));

Hopefully, I'm not bothering you with my questions regards.

[Star Strider]

My pleasure.

To do any statistics on your signal, such as autocorrelation, cross-correlation, entropy, and the rest, you likely need your entire time-domain signal. (I’m not certain what metrics you are using or how you are calculating them.)

I would not use the fft results, especially those plotted in figure(3), because that plotting technique irreversibly ‘loses’ about 75% of your data. (Half your signal is gone, and you cannot recover the real and imaginary parts of the half in the plot from the absolute value alone.)

[neamah al-naffakh]

thank you so much, but I am bit confused , because if I should use the time-domain signal to measure (energy, entropy and correlation),So, what is the benefit of converting the time-domain signal into the frequency domain signal?

I'm really sorry to keep asking but I just want to understand this point. I was thinking after converting the time domain into frequency domain signal ( then I use frequency domain signal as an input array in order to calculate energy, entropy, etc.

[Star Strider]

My pleasure.

That depends on how you calculate your statistics, especially entropy and correlation. I’ve used the time domain signal the few times I’ve needed to do that. There may be other techniques I’m not aware of.

The principal benefit of converting the time domain signal to a frequency domain signal is to understand the frequency components and to get the power spectral density. (The spectral energy in the signal is the fft as you have calculated and plotted it in figure(3). The power spectral density is the square of that.) Otherwise, knowing the frequency components of your signal is particularly important if you want to design filters.

[neamah al-naffakh]

Dear Star Strider , thank you so much for your patience, you are such a great person with my respect

[Star Strider]

As always, my pleasure.

[Star Strider]

neamah al-naffakh’s Comment (20:20 UTC) is duplicated here:

Dear Star Strider, please find the attached PDF file.

I have highlighted what the author did in order to convert his time domain signal to frequency- domain and then calculate the Amplitude in a specific frequency but I don't know how to do this in Matlab?

the array called Segments is the original data, could you help me to re-build the code please?

[Star Strider]

I'm moving our discussion back up here.

The paper appears fairly straightforward with respect to identifying the ‘features’. You first need to design a series of bandpass filters to match the frequencies in 4. on page 1082. There are several ways to design filters in MATLAB, a relativey easy one is designfilt. (I usually prefer to design my own, and my procedure is here: How to design a lowpass filter for ocean wave data in Matlab?)

I would take the RMS of each filter output and, as the paper describes in 3., concatenate them into a feature vector.

That is how I interpret the sections you highlighted.

I outlined one method of designing a filter bank in How to shape the spectrum of an audio waveform? so you can use it to get started.

Designing filters for your application requires 6 filters. If you attach a sample of your signal matrix as a .mat file, I can help you design your filters and test them. Your data must include the signals themselves and a time vector common to all of them.

I understand if you do not want to post your data online. If so, I only need to know the sampling frequency of your signal. I can help you design your filters, but would not be able to test them on your data. You would have to do that.

[neamah al-naffakh]

Dear Star Strider, I'm happy to share my data with you. Let me first explain to you what my data look like 1-the attached file (test2.csv) contains 4 columns, first col is the time, second col is the acceleration data of x-axes, third col is the acceleration data of y-axes, and the last col contains the the acceleration data of y-axes.

2- the first column in the file ( which represents the time), you will see a repeatation such as the number 0 repeats 32 or 31 times ( that means I got 32 samples of the movement data per second)

3- This is 5 min of acceleration data (which are collected at a rate of 30-32 samples per second for the x, y and z axes).

4- In order to remove noise, I have designed low pass filter for each axis with specific cut-off frequency

5- After the noise has been removed, the raw acceleration data has segmented into fixed windows, each window contains 10 seconds of the acceleration data ( so about 320 in each segment).

6- After data has been divided into segments, I have calculated some statistical features such as (mean, min, max, etc), but these features have extracted from the time domain signal.

7- All I need now is to convert the time domain into frequency domain and then extract extra features (Energy, Amplitude)

8- I have attached my code, please have a look and I will appreciate if you have any comments or correction on the code.

9- I have added comments to make the code clear and understandable.

10- please check your email, I have sent a special request and waiting your response if you dont mind.

Kind regards.

[Star Strider]

In order to do any meaningful signal processing on your data, it must all be sampled at the same constant and unvarying interval. If the sampling is 30-32 samples/second, it would be easiest initially to assign a time vector sampled at 1/31 second. It may not be necessary to resample it.

Ideally, your sampling time vector should have at least 3-digit precision. One-digit precision is simply not enough.

I didn’t see that the paper you quoted did any binning or segmentation of your data, but I didn’t look closely. You can certainly incorporate that if you want to.

My approach will first assign a constant-interval time vector, then to design the filters to extract the 6 ‘features’ mentioned. There is no need for any other pre-processing, since the filters themselves will remove high-frequency noise and d-c offset.

[neamah al-naffakh]

Hi Sir, i have updated my previous comment. for the sampling rate, as you can see the application sometimes output 33 samples per second, or 32,31,30 per second.

for the time interval, please see the attached file ( which is the original one without any modification)

about segmenting the data, it's my approach to segment the data into fixed segment size (10 seconds) and then extracting a set of time domain features as well as frequency domain features in order to use them for authentication regards.

[neamah al-naffakh]

this is an example of the segmentation approach ( called sliding window). you can see the highlighted part .

the only different between me and them I didn't use overlap window.

by the way, you are right we need to fix the sampling rate. but I'm wondering can we interpolate the sampling rate into 20, 25 or 30 per second?

2- in the attached paper, they have highlighted Bin ( which is N-Bin Histogram ),I don't know how to calculate it exactly in Matlab for each segment !

3- the most important features in this paper are MFCC, BFCC1, BFCC2. Could you please help me to extract them?

many regards and respect for your kind answer , I will never forget your help

[Star Strider]

Because your signal was sampled so infrequently (at about 31.8 samples/minute), it was necessary to resample them to a frequency high enough to permit the Nyquist frequency (the highest uniquely identifiable frequency in a sampled signal) to exceed the frequency of your highest filter.

This code reads your data, creates a unique time vector for the original signal, resamples the signal, designs the filters (a parallel bank of 6), filters your signal, and takes the RMS value of each signal in each band. It also plots your original signal, the resampled signal, the filter Bode plots (passbands), and the filtered signals. I did my best to document it internally with comments. (The code is not as efficient as it could be, because I left in intermediate steps to remind me what I did and why. It may also make it easier for you to understand what I did.)

The code:

[d,s,r] = xlsread('neamah al-naffakh test2.csv');
tx = d(:,1);                                                    % Original ‘Time’ Vector
smp_cts = accumarray(tx+1, 1);                                  % Samples/Minute Counts
Sm = std(smp_cts);
Fm = mean(smp_cts);                                             % Mean Sampling Frequency (samples/min)
Fs = Fm/60;                                                     % Mean Sampling Frequency (Hz)
Fn = Fs/2;                                                      % Nyquist Frequency
Ts = 1/Fs;                                                      % Sampling Interval (sec)

tv = linspace(0, size(tx,1), size(tx,1))/Fm;                    % ‘Synthetic’ Time Vector (min)
tvs = tv*60;                                                    % ‘Synthetic’ Time Vector (sec)

Fsd = 15;                                                       % Sampling Frequency Desired (Hz)
[N,D] = rat(Fsd/Fs);                                            % Integer Values For Resampling Rates

dr = resample([d(:,2:end)], N, D);                              % Resample Signals to ‘Fsd’ Hz
tvrm = linspace(0, tv(end), size(dr,1));                        % Resampled ‘Synthetic’ Time Vector (min)
tvr = tvrm*60;                                                  % Resampled ‘Synthetic’ Time Vector (sec)
Tsr = mean(diff(tvr));                                          % Resampled Signal Sampling Interval
Fsr = 1/Tsr;                                                    % Resampled Signal Sampling Frequency
Fnr = Fsr/2;                                                    % Resampled Signal Nyauist Frequency

figure(1)
plot(tv', d(:,2:end))
grid
axis([0  10    -1  2])
title('Original Signals')

figure(2)
plot(tvrm', dr)
grid
axis([0  10    -1  2])
title('Resampled Signals')

BandsMtx = [0.5:5.5;  1.5:6.5]';                                % Filter Passband Matrix

Rp =  1;                                                        % Passband Ripple
Rs = 40;                                                        % Stopband Ripple

figure(3)
for k1 = 1:size(BandsMtx,1)                                     % Design Filters, Display Bode Plots
    Wp = BandsMtx(k1,:)/Fnr;                                    % Normalised Passband
    Ws = (BandsMtx(k1,:) + [-0.2  0.2])/Fnr;                    % Normalised Stopband
    [n,Ws] = cheb2ord(Wp, Ws, Rp, Rs);                          % Use Chebyshev Type II Bandpass Filters
    [b,a] = cheby2(n,Rs,Ws);
    [sos{k1},g(k1)] = tf2sos(b,a);                              % Second-Order-Section For Stability
    freqz(sos{k1}, 4096, Fsr)
    if k1 == 1
        hold on
    end
end
hold off

for k1 = 1:size(BandsMtx,1)
    Band{k1} = filtfilt(sos{k1},g(k1),dr);                      % Filter Signals
    RMS_Band(k1,:) = sqrt(mean(Band{k1}.^2));                   % Calculate Signal RMS Values
end

figure(4)
for k1 = 1:size(BandsMtx,1)
    subplot(6,1,k1)
    plot(tvrm, Band{k1})
    grid
    axis([0  10    ylim*0.4])
    title(sprintf('Band %.0f: [%.1f %.1f] Hz',k1, BandsMtx(k1,:)))
end

If you could get your time vector in milliseconds instead of minutes, it would make the analysis much easier.

[neamah al-naffakh]

Thank you so much Sir, i will read it and come back again if you dont mind, please let me know if you want to end the discussion because your comments really help me.

Sir, have you read my previous comment please? about MFCC, BFCC1, and BFCC2,because to be honest when it comes to frequency domain features i dont really understand that much

could you help me to extract MFCC, BFCC1, and BFCC2 based on the paper in this comment ?

[Star Strider]

My pleasure.

I won’t end it until I can no longer offer you substantive help. There are limits to my knowledge, expertise, and time, so providing we remain within those constraints, I’ll help as I can.

Please have a go with my code. It took me a few hours to design the time vectors for your original data, then resample them and design, implement, and test the filters on them. My code works, and should provide what you need.

I’ll explore the paper with MFCC, BFCC1, and BFCC2 later, and let you know if I can help you with them.

--------------------------------------------------------------

EDIT — 27 Aug 2016 14:43 UTC

The mel-cepstral and bark-cepstral analyses are areas I would prefer not to explore just yet. However there are online resources referenced by [11.] in the paper by Nickel, et al., particularly Reproducing the feature outputs of common programs using Matlab and melfcc.m and PLP and RASTA (and MFCC, and inversion) in Matlab using melfcc.m and invmelfcc.m. I encourage you to explore them. Since the code is already written, all you need to do is to understand how they work and how to use them. (I am bookmarking them for my later reference, to explore when I have time.)

[neamah al-naffakh]

Hi Sir,please do not upset from this comment :(.

Although you have documented your code but I don't understand it!! and don't understand what are you trying to do?

let me explain some points

1- the first column in (test2.csv) is the time vector in seconds not minute ( e.g. all the records have recorded in one second, particularly, second 16) see this

2- what are trying to do on this code and why?

    % smp_cts = accumarray(tx+1,1);                                  % Samples/Minute Counts
Sm = std(smp_cts);
Fm = mean(smp_cts);                                             % Mean Sampling Frequency (samples/min)
Fs = Fm/60;                                                     % Mean Sampling Frequency (Hz)
Fn = Fs/2;                                                      % Nyquist Frequency
Ts = 1/Fs;                                                      % Sampling Interval (sec)

3- again what are you trying to do here?

   tv = linspace(0, size(tx,1), size(tx,1))/Fm;                    % ‘Synthetic’ Time Vector (min)
tvs = tv*60;                                                    % ‘Synthetic’ Time Vector (sec)
Fsd = 15;                                                       % Sampling Frequency Desired (Hz)
[N,D] = rat(Fsd/Fs);                                            % Integer Values For Resampling Rates

    dr = resample([d(:,2:end)], N, D);                              % Resample Signals to ‘Fsd’ Hz
tvrm = linspace(0, tv(end), size(dr,1));                        % Resampled ‘Synthetic’ Time Vector (min)
tvr = tvrm*60;                                                  % Resampled ‘Synthetic’ Time Vector (sec)
Tsr = mean(diff(tvr));                                          % Resampled Signal Sampling Interval
Fsr = 1/Tsr;                                                    % Resampled Signal Sampling Frequency
Fnr = Fsr/2;

A- First, I want to design a filter to remove the noise from the original data and store the filtered signal in an array.

B- Because the data has recorded in different sampling rate (e.g. sometimes I get 32,33,31,30 samples or values per second, I want to get fix sample rate such as 30 samples per second).

let's do this first and then i will go to the next question

let me know if you are still happy to keep the discussion

with my respect

[Star Strider]

1. Your data are sampled at about 31 times per minute, or about 0.6 Hz. I resampled it to 15 Hz, because the upper limit on your filters is 6.5 Hz. That means you need a Nyquist frequency of at least 13 Hz.

2. The accumarray call is an easy way to see how regularly your data are sampled. For example the first minute has 32 samples, and some others have 31. This gives me an idea of what the time vector should be, so that I can reconstruct it before resampling it.

3. That section of the code creates a correct time vector (in seconds) for your original data, resamples your data, then recalculates a second time vector for the resampled data. The plots (in figure(1) and figure(2)) demonstrate that all this is done correctly.

4. Describing in detail:

A. A separate filter to remove high-frequency noise and baseline drift is not necessary here. The first paper you provided divides the signal into 6 bands, each bandpass filter removing baseline drift and high-frequency noise individually.

This section of code designs those filters and plots their characteristics:

BandsMtx = [0.5:5.5;  1.5:6.5]';                                % Filter Passband Matrix

Rp =  1;                                                        % Passband Ripple
Rs = 40;                                                        % Stopband Ripple

figure(3)
for k1 = 1:size(BandsMtx,1)                                     % Design Filters, Display Bode Plots
    Wp = BandsMtx(k1,:)/Fnr;                                    % Normalised Passband
    Ws = (BandsMtx(k1,:) + [-0.2  0.2])/Fnr;                    % Normalised Stopband
    [n,Ws] = cheb2ord(Wp, Ws, Rp, Rs);                          % Use Chebyshev Type II Bandpass Filters
    [b,a] = cheby2(n,Rs,Ws);
    [sos{k1},g(k1)] = tf2sos(b,a);                              % Second-Order-Section For Stability
    freqz(sos{k1}, 4096, Fsr)
    if k1 == 1
        hold on
    end
end
hold off

This section of code filters your 3 signals into the six bands and stores them in the cell array ‘Band’, one cell for each band, calculates the RMS value for each band and stores it in ‘RMS_Band’:

for k1 = 1:size(BandsMtx,1)
    Band{k1} = filtfilt(sos{k1},g(k1),dr);                      % Filter Signals
    RMS_Band(k1,:) = sqrt(mean(Band{k1}.^2));                   % Calculate Signal RMS Values
end

If you instead want to do only one filter, decide what frequencies you want to include as your signal, and use the filter design procedure in my loop (but without the loop since you’re only designing one filter) to specify the passbands. (I did not see any noise in your signal, only what I would consider to be normal movement.) Filter my resampled signal, not your original signal, because your original signal is undersampled. Note that the upper passband of your bandpass filter is limited to at most 7.4 Hz.

B. I already did that in the original and resampled signals.

C.

D. You can do that on my resampled data. My data are resampled to a frequency of 15 Hz (samples/second).

E. You can do that on my resampled data.

F. The other paper you posted suggested using a Support Vector Machine for the classification. (I also found and posted a number of sources for MATLAB routines to do the mel-cepstral analyses.) I’ve classified electroencephalogram data but not recently. Classification is quite involved, and something I would prefer not to get into.

[neamah al-naffakh]

Hi Sir, sorry to bother you again, Sir, it seems that I made you confused .

please see the attached file ( the first column, it the Date Time which has the form of Hour, min, second ( start with this 10:34:16 PM ) 10-- hour , 34 --min, 16 --second

you can see that the data has recorded at a sampling rate between 33-30 per second not min !

now, forget about the paper that I have attached let me explain to you step by step and it is up to you Sir if you don't want to help because I know that I bothered you so much

1- the first preprocessing step to my original data, I want to linearly interpolate the data in order to get a fixed sampling rate(e.g. 30 samples per seconds) and if possible, could you write a comment in the code to help me if I want to change the sampling rate into 20 or 40 samples per second? because I want to test more than one setting for the sampling rate and select the best.

2- After interpolation, I want to design a filter to remove the noise from the original data and store the filtered signal in an array ( please make sure the array ,which contains the filtered data, has four columns ( the first column contains time in second , the second, third and fourth columns contain the data of x,y,z-axes respectively)

could we do these two steps first and ignore the paper that i have sent please?

many regards and sorry for the inconvenience that I caused for you

with my respect

[Star Strider]

The first record you attached (as test2.csv) did not have any distinct time vector at all. I had no idea what the time base was, and I assumed minutes.

There are 32 identical values of 10:34:16 PM. That means that the first second of your signal is sampled 32 times. Your entire signal has an average sampling frequency of about 31.85 Hz, or one sample every 0.031 seconds.

I resampled your signal at 45 Hz, since that makes the filter design easier. (I also tested it at 30 Hz, and several other frequencies, 15 Hz and above. It ran with all without error.) You can specify:

Fsd = 30;

and the code will run and give the correct results. Note that because of the filter passbands, 15 Hz is the lowest usable sampling frequency.

The code (otherwise essentially unchanged) with the redefined time base:

[d,s,r] = xlsread('neamah al-naffakh test2.csv');
tx = d(:,1);                                                    % Original ‘Time’ Vector
smp_cts = accumarray(tx+1, 1);                                  % Samples/Second Counts
Sm = std(smp_cts);
Fm = mean(smp_cts);                                             % Mean Sampling Frequency (samples/sec)
Fs = Fm;                                                        % Mean Sampling Frequency (Hz)
Fn = Fs/2;                                                      % Nyquist Frequency
Ts = 1/Fs;                                                      % Sampling Interval (sec)

tv = linspace(0, size(tx,1), size(tx,1))/Fm;                    % ‘Synthetic’ Time Vector (sec)

Fsd = 45;                                                       % Sampling Frequency Desired (Hz)
[N,D] = rat(Fsd/Fs);                                            % Integer Values For Resampling Rates

dr = resample([d(:,2:end)], N, D);                              % Resample Signals to ‘Fsd’ Hz
tvr = linspace(0, tv(end), size(dr,1));                         % Resampled ‘Synthetic’ Time Vector (sec)
Tsr = mean(diff(tvr));                                          % Resampled Signal Sampling Interval (sec)
Fsr = 1/Tsr;                                                    % Resampled Signal Sampling Frequency (Hz)
Fnr = Fsr/2;                                                    % Resampled Signal Nyauist Frequency (Hz)

figure(1)
plot(tv', d(:,2:end))
grid
axis([0  10    -1  2])
title('Original Signals')

figure(2)
plot(tvr', dr)
grid
axis([0  10    -1  2])
title('Resampled Signals')

BandsMtx = [0.5:5.5;  1.5:6.5]';                                % Filter Passband Matrix

Rp = 10;                                                        % Passband Ripple
Rs = 40;                                                        % Stopband Ripple

figure(3)
for k1 = 1:size(BandsMtx,1)                                     % Design Filters, Display Bode Plots
    Wp = BandsMtx(k1,:)/Fnr;                                    % Normalised Passband
    Ws = (BandsMtx(k1,:) + [-0.3  0.3])/Fnr;                    % Normalised Stopband
    [n,Ws] = cheb2ord(Wp, Ws, Rp, Rs);                          % Use Chebyshev Type II Bandpass Filters
    [b,a] = cheby2(n,Rs,Ws);
    [sos{k1},g(k1)] = tf2sos(b,a);                              % Second-Order-Section For Stability
    freqz(sos{k1}, 4096, Fsr)
    if k1 == 1
        hold on
    end
end
hold off

for k1 = 1:size(BandsMtx,1)
    Band{k1} = filtfilt(sos{k1},g(k1),dr);                      % Filter Signals
    RMS_Band(k1,:) = sqrt(mean(Band{k1}.^2));                   % Calculate Signal RMS Values
end

figure(4)
for k1 = 1:size(BandsMtx,1)
    subplot(6,1,k1)
    plot(tvr, Band{k1})
    grid
    axis([0  10    -0.5  0.5])
    title(sprintf('Band %.0f: [%.1f %.1f] Hz',k1, BandsMtx(k1,:)))
end

[neamah al-naffakh]

Hi sir, thank you so much.

I am going through your code now, I understand what are you doing in the first part of your code.

I guess it is easy to you and me as well, if I ask one question at the time to prevent the misunderstanding,

so I'm wondering is the above code to re-sample my data at a fixed rate and design a filter to remove the noise from the signal? or do something else?

because i want to forget what the author did in his paper, i want to follow my approach and you help me so much in re-sample my data.

regards.

Many respect

[Star Strider]

My pleasure.

You need to resample your signal to a higher sampling frequency if you want to use a filter with a passband or stopband frequency above 15 Hz. If you use my code, you do not need to use a separate filter to remove the baseline drift or baseline offset or high-frequency noise. The ‘Band’ filters do that themselves.

If you want to create a separate filter, you can use the filter design procedure I used to design the band filters in the loop. Use the fft of your signal to decide on the passbands and stopbands. To design the filter, you would just change the code to use these assignments:

[n,Ws] = cheb2ord(Wp, Ws, Rp, Rs);                          % Use Chebyshev Type II Bandpass Filters
[b,a] = cheby2(n,Rs,Ws);
[sos,g] = tf2sos(b,a);                                      % Second-Order-Section For Stability 
freqz(sos, 4096, Fs)                                        % ‘Fs’ Is The Sampling Frequency Of Your Signal

and then to filter the signal:

FilteredSignal = filtfilt(sos,g,Signal);                    % Filter Signals

You have to decide what the passband, Wp, and stopband, Ws frequencies are. I would start by using my values for Rp and Rs, although you may want to experiment with them.

[neamah al-naffakh]

Sir,you are really so kind!! I have never seen someone patience like you.

I will go through your code and comments and let you know if I have a question. please keep in touch with me please.

really appreciate your help.

[Star Strider]

My pleasure.

[neamah al-naffakh]

Hi Sir, I hope you are ok .

As you know that I have divided the row data into segments, and each segment contains 300 values now.

I would like to calculate 10-Bin-Histogram for each segment, some of them call it Binned Distribution.

simple description about it

determining the range of values for each axis (maximum – minimum), divide this range into 10 equal sized bins, and then record the fraction of the 300 values that fall within each of the bins.

I really appreciate your kind help!

I did that but it seems not correct

    hist_arr_x{nn,1}=histcounts(X_Segments,10);
     hist_arr_y{nn,1}=histcounts(Y_Segments,10);
     hist_arr_z{nn,1}=histcounts(Z_Segments,10);

Note

X_Segments, contains 300 values of X

Y_Segments ,contains 300 values of Y and

Z_Segments contains 300 values of Z

[Star Strider]

The code looks correct to me. If you want to know the edges of the bins, you have to request them:

[hist_arr_x{nn},edges_x{nn}] = histcounts(X_Segments,10);

[Star Strider]

‘but does that code meet the description of 10-Bin-histogram’

It does.

Example:

x = randi(99, 1, 300);                                                  % Create Data
[cts,edg] = histcounts(x, 10);                                          % Histogram Counts & Edges
cts_size = numel(cts)                                                   % Confirm Number Of Bins

bin_ctrs = edg(1:end-1) + mean(diff(edg))/2                             % Calculate Bin Centres (For Plot)

cts_size =
    10

bin_ctrs =
     5    15    25    35    45    55    65    75    85    95

[neamah al-naffakh]

ok Sir,

now, I got what you mean, but each segment has a different edge! for example,*segment1*(min=0, max=1.4), segment2(min=0.7, max=1.3). I want to fix the range of all segments(min and max).

and then calculate the 10-Bin-Histogram for each segment.

[Star Strider]

You can specify the edges. See the documentation for Specify Bin Edges.

Remember that for 10 bins, you need to specify 11 bin edges. They do not have to be regularly-spaced, but if they are not, my code for calculating the bin centres (the ‘bin_ctrs’ assignment) will not be correct.

[neamah al-naffakh]

Hi sir, i have calculated the (Min-Max) edges for all segments, and they are in the range between 0.56 to 1.41 (the Min Edge is 0.56 and the Max Edge is 1.41 )

how can specify Specify 10-Bin Edges.?

[Star Strider]

You specify 11 bin edges to get 10 bins with histcounts. The easiest way is to use the linspace function:

bin_edges = linspace(0.56, 1.41, 11);

[neamah al-naffakh]

Hi sir, thank you so much. i have done what you advised.

do you mind to open another topic and keep the discussion? as i have another questions related to my data

many regards and thanks for your kind help

[Star Strider]

My pleasure.

It depends what the topic is.

Consider asking a new Question in a new thread if it is actually a ‘new topic’.

[neamah al-naffakh]

Hi sir, thank you again.

it's same but different question

because now I have done the calculation of some features from the time domain signal, and I want to derive some features from the frequency domain signal.

Therefore, I have written new topic to make it easily readable and useful for reader.

I will keep this one for time domain signal only

here is the new topic

[neamah al-naffakh]

Dear sir, thanks for your support

I have fixed the time stamp in my application( so the data now is collected in milliseconds). hence, that means I have distinct time vector or time stamp in milliseconds (see the attached file).

I am wondering if we can resample my data at 30 Hz by using your previous code? because when I run your code I get this error.

*Error using accumarray Requested 1460000000001x1 (10877.8GB) array exceeds maximum array size preference.

Creation of arrays greater than this limit may take a long time and cause MATLAB to become unresponsive. See array size limit or preference panel for more information.

Error in AccGyr (line 4) smp_cts = accumarray(tx+1, 1);*

what can we change in this code in order to make it run correctly with the attached file ?

    % [d,s,r] = xlsread('test2.csv');                                 % (d) contains the original data 
tx = d(:,1);                                                    % Original ‘Time’ Vector
smp_cts = accumarray(tx+1, 1);                                  % Samples/Second Counts
Sm = std(smp_cts);
Fm = mean(smp_cts);                                             % Mean Sampling Frequency (samples/sec)
Fs = Fm;                                                        % Mean Sampling Frequency (Hz)
Fn = Fs/2;                                                      % Nyquist Frequency
Ts = 1/Fs;                                                      % Sampling Interval (sec) means how many seconds is taken to record one sample
tv = linspace(0, size(tx,1), size(tx,1))/Fm;                    % ‘Synthetic’ Time Vector (sec) (creat new {distinct} time vector)

    % Fsd = 30;                                                       % Sampling Frequency Desired (Hz)
[N,D] = rat(Fsd/Fs); 
dr = resample([d(:,2:end)], N, D);                              % (dr) (data_re-sampled)contains the data after it Resample to the New ‘Fsd’ Hz ( the new Sampling Rate)
tvr = linspace(0, tv(end), size(dr,1));                         % Resampled ‘Synthetic’ Time Vector (sec)
Tsr = mean(diff(tvr));                                          % Resampled Signal Sampling Interval (sec)
Fsr = 1/Tsr;                                                    % Resampled Signal Sampling Frequency (Hz) make the sample rate 30 
Fnr = Fsr/2;     

    % tvr=tvr';
sampled_arr(:,1)=tvr' ;         % tvr (time_vector_re-sampled) the new time vector
sampled_arr(:,3:5)=dr;          % dr (data_re-sampled) new resample data,

[Star Strider]

My pleasure.

No file was attached. However if you have a time vector in milliseconds, it will probably be relatively straightforward to resample it to 30 Hz. I would have to see your data.

[neamah al-naffakh]

hi sir, sorry i forgot

by the way, i can't attach the file because it's 21MB so i have uploaded here

[Star Strider]

It might be easier if you just uploaded the first 2 MB here. I just need enough to test with.

[neamah al-naffakh]

hi sir, please see the attached file regards.

[Star Strider]

We still have a problem. The timestamps are not unique. The first 5 rows of the data in your Excel file are:

data =

    Columns 1 through 3

               1470000000000                       NaN              -0.022460938
               1470000000000                       NaN              -0.013183594
               1470000000000                       NaN              -0.016601563
               1470000000000                       NaN              -0.006835938
               1470000000000                       NaN               -0.05810547

    Columns 4 through 5

                  -0.7216797                -0.7331543
                  -0.6218262               -0.72021484
                  -0.6052246                 -0.685791
                  -0.5644531                -0.7246094
                  -0.5559082               -0.68408203

The 'DateTime' strings are not unique either. I checked with the original Excel file and there were no changes and no truncation of any of the numbers in importing them.

There seem to have not been any significant changes in the data.

[neamah al-naffakh]

hi sir,

The timestamps are unique but it seems when i uploaded only 5 MB from the file the time stamp has been changed

could you please download this file

[neamah al-naffakh]

Hi sir,

here we go please see the attached file.

regards.

[Star Strider]

Great!

The times are definitely unique now. When I print the first 5 rows, I get:

time =
     1.473948671396000e+12
     1.473948671464000e+12
     1.473948671478000e+12
     1.473948671529000e+12
     1.473948671546000e+12

The first column is in milliseconds, correct? So we have to multiply the first column by 1E-15 to get seconds. When I do that, I get a sampling frequency of 30E+12 Hz, which seems a bit high.

I need to understand the units your time is in before I go further.

[neamah al-naffakh]

hi sir thanks for your kind help,

actully, I dont know the time units because this application is written by developers, I will back later to discuss this point with you

regards.

[Star Strider]

My pleasure.

We really need to know the time units. There are several different time systems in computers, so knowing that will help as well. For example, when did the time reference begin (when did the time = 0)? A datetime object might make the time conversion trivial, or we might have to write code to do the conversion.

Answer 2

[Marcel]

Hint: In matlab plot(fft(y)) will plot only the real part of your spectrum there by ignoring the imaginary values in each of the frequency samples. You may do the DFT naively by plotting the amplitude of your spectrum or the phase i.e. plot(abs(fft(y))) etc. You may want to look for the pwelch or periodogram accordingly to what you want (PSD ....).

[neamah al-naffakh]

Hi Marcel, thank you so much for your answer as well I will consider your answer , but could you please keep in touch with me in order if I have a question to your comment?

[Marcel]

Cool no problem.

[neamah al-naffakh]

Hi Marcel, I tried to calculate the Magnitude , Phase, and Energy for the frequency domain as Magnitude can be calculated

XR Is the real part while xr is imaginary part signal.

so,I have used this code,

1- FD_Signal_XYZ=fft(Segments);  % To convert the time domain into frequency domain 
2- magnitude_XYZ =abs(FD_Signal_XYZ);  % To measure the Ampltiude in the frequency signal      
3-phase_XYZ=unwrap(angle(FD_Signal_XYZ{1}));        
4-Energy_XYZ = (FD_Signal_XYZ).^2;
5-total_Energy_XYZ = sum(Energy_XYZ);

howeover, the Output of the magnitude and Phase were many values? I'd like to get one value for both of them for the selected freqency signal.

Many regards.

[Star Strider]

Every frequency in a Fourier-transformed signal has a single magnitude and phase.

The ‘energy’ at a particular frequency is the amplitude of the signal at that frequency. This line in your code calculates the power (amplitude squared) at every frequency:

Energy_XYZ = (FD_Signal_XYZ).^2;

If you want the direct-current equivalent of the energy in your entire signal, take the root-mean-squared value of the entire signal:

RMS_Amplitude = sqrt(mean(Segments.^2));

The equivalent d-c power is the square of that:

RMS_Power = RMS_Amplitude^2;

[neamah al-naffakh]

dear Star Strider, you are always helpful and so kind with me thank you so much.

I'm sorry to ask again but let me ask you some questions related to your answer

1-do you mean I should calculate the energy for the time domain signal( which is stored in the array called Segments not for output of FFT function)?

if yes, actually, I already calculted root-mean-squared for the time domain signal but I am looking to calculate some measurement for the frequency domain signal

2- what do you mean by d-c power?

[Star Strider]

My pleasure.

1. Yes. It seems you have already calculated the energy and power of each frequency in your Fourier-transformed frequency domain signal.

2. If you rectify and filter an alternating-current (a-c) signal, the equivalent mean energy is the RMS value of the a-c signal. This is the direct current (d-c) equivalent energy. The d-c power is the square of the d-c energy.

At least, that is how I learned it.

[neamah al-naffakh]

Dear Star Strider, do you mind if i attach a PDF file that illustrates what I want to do with the frequency signal?

as the the paper exactly explained what I'd like to measure, please let me know if you are happy with that

With My Respect

[neamah al-naffakh]

Dear Star Strider, please find the attached PDF file.

I have highlighted what the author did in order to convert his time domain signal to frequency- domain and then calculate the Amplitude in a specific frequency but I don't know how to do this in Matlab?

the array called Segments is the original data, could you help me to re-build the code please?

Answer 3

[RAUNAK GUPTA]

please help me I want to design audio control mood lamp.....in which I have to record an real time audio data in matlab and interfacing it with microcontroller I have to control light.....HOW to convert a real time audio data in frequency domai?please help me with the code

matlab code

clear all; clc;
rec=audiorecorder(8000,8,1);
ard = arduino('com4');
a = input('Press 1 to start recording\n');
recordblocking(rec,20);
data=getaudiodata(rec);
play(rec);
plot(data);
Fs=8000;
N = length(data);
datafft=fft(data);
datafft_abs=abs(datafft/N);
datafft_abs=datafft_abs(1:N/2+1);
f=Fs*(0:N/2)/N;
arduin(ard,f);
figure;
plot(f,datafft_abs);
clearvars ard;

arduin funtion

function [] = arduin(a,freq)
 %This is self made function.
 %This function is made for using arduino.
 %This function lits up the LEDs according to frequencies it get.

fr = 2; %Frequnecy to be set for red light
fb = 3; %Frequnecy to be set for blue light
configurePin(a,'D8');
configurePin(a,'D9');
configurePin(a,'D13');
for i=1:length(freq)
    f = freq(i);
    if (f < fr)
        disp('RED');
        writeDigitalPin(a, 'D9', 0);
        writeDigitalPin(a, 'D13', 0);
        writeDigitalPin(a, 'D8', 1);
    else
        if (f<=fb) && (f=fr)
            disp('GREEN');
            writeDigitalPin(a, 'D8', 0);
            writeDigitalPin(a, 'D13', 0);
            writeDigitalPin(a, 'D9', 1);
        else
            if (f>=fb)
                disp('Blue');
                writeDigitalPin(a, 'D9', 0);
                writeDigitalPin(a, 'D8', 0);
                writeDigitalPin(a, 'D13', 1);
            else
                disp('Frequency out of range');
                writeDigitalPin(a, 'D8', 1);
                writeDigitalPin(a, 'D9', 1);
                writeDigitalPin(a, 'D13', 1);
            end
        end
    end
end