**You are now following this question**

- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.

# convert the time domain signal into frequency domain signal

41 views (last 30 days)

Show older comments

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

### Accepted Answer

Star Strider
on 18 Aug 2016

##### 52 Comments

neamah al-naffakh
on 18 Aug 2016

Star Strider
on 18 Aug 2016

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
on 19 Aug 2016

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
on 19 Aug 2016

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
on 21 Aug 2016

Edited: neamah al-naffakh
on 21 Aug 2016

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
on 21 Aug 2016

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
on 21 Aug 2016

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
on 21 Aug 2016

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
on 21 Aug 2016

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
on 21 Aug 2016

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
on 21 Aug 2016

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
on 21 Aug 2016

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
on 24 Aug 2016

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

Star Strider
on 25 Aug 2016

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
on 25 Aug 2016

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
on 26 Aug 2016

Edited: neamah al-naffakh
on 26 Aug 2016

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
on 26 Aug 2016

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
on 26 Aug 2016

Edited: neamah al-naffakh
on 26 Aug 2016

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
on 26 Aug 2016

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
on 26 Aug 2016

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
on 26 Aug 2016

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

Star Strider
on 26 Aug 2016

Edited: Star Strider
on 27 Aug 2016

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
on 27 Aug 2016

Edited: neamah al-naffakh
on 10 Sep 2016

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
on 27 Aug 2016

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
on 28 Aug 2016

Edited: neamah al-naffakh
on 28 Aug 2016

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
on 28 Aug 2016

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
on 28 Aug 2016

Edited: neamah al-naffakh
on 28 Aug 2016

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
on 28 Aug 2016

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
on 28 Aug 2016

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.

neamah al-naffakh
on 4 Sep 2016

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
on 4 Sep 2016

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
on 4 Sep 2016

‘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
on 4 Sep 2016

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
on 4 Sep 2016

neamah al-naffakh
on 5 Sep 2016

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
on 5 Sep 2016

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
on 5 Sep 2016

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
on 5 Sep 2016

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
on 6 Sep 2016

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

neamah al-naffakh
on 14 Sep 2016

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
on 14 Sep 2016

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.

Star Strider
on 15 Sep 2016

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

neamah al-naffakh
on 15 Sep 2016

Star Strider
on 15 Sep 2016

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
on 15 Sep 2016

hi sir,

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

neamah al-naffakh
on 15 Sep 2016

Edited: neamah al-naffakh
on 15 Sep 2016

Hi sir,

here we go please see the attached file.

regards.

Star Strider
on 15 Sep 2016

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
on 16 Sep 2016

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
on 16 Sep 2016

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.

### More Answers (2)

Marcel
on 24 Aug 2016

##### 8 Comments

neamah al-naffakh
on 24 Aug 2016

neamah al-naffakh
on 25 Aug 2016

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
on 25 Aug 2016

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
on 25 Aug 2016

Edited: neamah al-naffakh
on 25 Aug 2016

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
on 25 Aug 2016

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
on 25 Aug 2016

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
on 25 Aug 2016

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?

RAUNAK GUPTA
on 15 Oct 2017

Edited: Walter Roberson
on 21 Nov 2021

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

### See Also

### Community Treasure Hunt

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

Start Hunting!**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list

How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

Americas

- América Latina (Español)
- Canada (English)
- United States (English)

Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)

Asia Pacific

- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)