texts: How can I detect potholes using an accelerometer mounted in a vehicle?

texts:
I am currently studying DSP and FFT, I am very new to this and have been doing electronics for a long time with Arduino and hobby projects. Recently, I am doing a project with the goal to map potholes along your daily commute in the car, measuring road roughness. This is not the road profile, but the subjected roughness of the ride the driver feels during the trip. I have an acceleration mounted to measure Z-axis "vertical" acceleration of the car, when a pothole is struck the shocks and springs dampen the force as per the quarter car model for simplicity.

Basically I want to be able to create a detector that will detect potholes by using FFT to find the dominating frequency of the vertical acceleration, then I can maybe match it to pattern to detect door slamming of the car door, or a man-hole cover from a pothole maybe in the future.

I'm not sure if FFT is the way to go, any input or advice and ideas are greatly appreciated, I have made a low pass filter to filter out the higher frequency vibration from the drive train and engine.

Any ideas on the subjected "roughness" felt by the passengers/driver in the car? I was thinking of doing road tests on a smooth road, gravel road, brick road, bad road, and coming up with a threshold of the magnitude of the vertical acceleration and making a scale somehow on what is considered "rough"

The potholes detected would be also logged with GPS their locations, I'm trying to use the Rasberry Pi for this, with data storage and its speed as a good project for it.

Thanks any help or ideas is appreciated, books, tutorials, wisdom, etc

UPDATE:

Carry over additional data here:

https://electronics.stackexchange.com/questions/56238/accelerometer-data-smoothing-filtering-pothole-detection


texts: What is the reasoning behind the deviration of propogation of uncertainty?

texts: When considering the uncertainty of a signal which is determined by multiple inputs propagation of uncertainty states that for measurement $$ y = f(x), x=\{x_1, x_2,..., x_N\} $$ uncertainty in $y$ is $$u_y=\sum_i \frac{\partial f}{\partial x_i} u_{x_i} $$ My question is how do we get to this rule. I know that it can be obtained by a taylor expansion of $f_i$ but I don't see why this is a sensible/logical thing to do.


texts: Regarding Bode plots; $H(s)$ and $H(j\omega)$

texts:
In circuit analysis, I understand the use of Laplace Transforms to obtain the impedance of a linear RLC circuit, ie transforming from the time domain to the frequency domain. In most texts I have seen and classes I have taken, the Bode plot looks at the rational transfer function of $H(j\omega)$, which is effectively the Fourier Transform, and not $H(s)$.

Since $\omega$ is considered to be the variable that represents frequency, what is $\sigma$ if $s = \sigma + j\omega$? Does $\sigma$ have a physical interpretation?


texts: Normalize two or more sound signals in order to let them with the same duration time and amplitude

texts:
I am very new in signal processing area and I am trying to solve the task in order to make an experiment with audio files from different languages, so, I am trying to let some ogg audio files with the same amplitude(wave height) and playing duration, so that the first wave can have the same amplitude from the second one.duration
My file:

Original file(The file that I am wanting to convert the previous file to the format of this file):

Any ideas in how could I begin on it? I am using Octave, or any algorithm to it in any other language (Python, MATLAB, Scala)?


texts:
I wish to test whether two time series are equal. So I believe best way to define equivalence is that given two time series, say $\{x1_t\}$ and $\{x2_t\}$, we show that both the series come from the same stochastic process, $\{X_t\}$.

How to do this?


texts:
Recently, I've read several papers where zero-lag Butterworth filters are used (both low-pass and high-pass).

I have a fair understanding of (Butterworth) filters, but how can we design one with zero lag? Can someone explain me how this works?

From my understanding, any low-pass filter yields an intrinsic delay.

Or does this involve bidirectional filtering?


texts: Whittaker-Shannon ($\mathrm{sinc}$) interpolation for a finite number of samples

texts:
Given an  infinite number of samples $(N)$, a higher (or lower) number of samples $(cN)$ can be derived using sinc interpolation followed by sampling. How can this be applied to finite length signals?

With $\mathrm{sinc}$ interpolation, one can derive a continuous-time signal as:

$$y(t) = \sum^{\infty}_{n=-\infty} y[n]\mathrm{sinc}\left({t\over T}-n\right)$$

For a finite number of sample points, should (can) we consider the $x[n]$ in the picture as

$$y[n] =
\begin{cases} x[n],  & \text{if } n \in [0, N-1] \\0, & \text{otherwise} \end{cases} $$

Or should $y[n]$ be considered as a periodic version of $x[n]$? (This link briefly addresses the same. The original stated form cannot be directly used with periodic signals)

$$y[n] = x\left[n\pmod N\right]$$

In the first consideration, outside the region $[0,\ N-1]$, if I understand correctly, the Gibb's phenomenon would result in a ringing effect. Would this completely invalidate any values predicted outside the non-zero region or is it only that the degree of inconsistency is high? (More specifically for points close to but just outside the boundary in the interpolated continuous-time signal)

I am interested to know whether the addition of zeros would pollute the input set of points during the interpolation stage.


texts: what should be a good bibliography for deep learning in 2020

texts: What should be a good bibliography for having a good overview of deep learning (on image/signal processing) in 2020 ?


texts: I could infer from books and wiki that a matching pursuit finds signal approximations from a redundant dictionary. Let f be the signal, ai be the basis and Ri be the residue.. here we choose the inner product of f and ai to be max such that Ri is less. What do we do next. what do we iterate?


texts:
I would like to create an algorithm that, given an iddata,  estimates different armax models with different na, nb, nc, nk, and finally show the compare plot.

I know how to use the armax estimation function, but how can I store the model in a variable in a such that I can compare them later?


texts: Audio distortion while implementing a low-pass filter

texts:
I attempted to implement a real-time audio filtering in a callback function using the simplest moving-average filter with coefficients [0.5 0.5] :

void process(void* in_, void* out_, DWORD length, void* user) {

    //data is float, mono
    int n = length/sizeof(float);
    float* in = (float*)in_;
    float* out = (float*)out_;
    float* buff = (float*)user;

    //buff holds the last input sample from the previous block
    //if this is a first block, it is 0
    out[0] = 0.5*in[0] + 0.5*buff;

    for (int i = 1; i < n; ++i)
    {
         out[i] = 0.5*in[i] + 0.5*in[i - 1];
    }

    buff[0] = in[n-1];
}

however, i get a distortion in the output, like an aliasing of some sort or clipping. But, to my knowledge none of this is possible to occur in this case. I've implemented the same code in Matlab ( although not in real-time ) and everything was fine.

I am using a BASS library.

I would like to know if i am missing something obvious here, cause i am running out of ideas. Thanks in advance!

EDIT: Thanks everyone for the answers. The above function is in the form for simplicity, but in actual from it looks like below (there is no output buffer, filtering must be done in-place) but it does exact same thing, although I will write it for convenience:

void CALLBACK filterDSPMono(HDSP handle, DWORD channel, void* buffer, DWORD length, void* user)
{
    float* x = (float*) buffer;
    int n = length / sizeof(float);
    float* buff = (float*) user;
    float temp;

    int i = 0;
    for (; i < n; ++i)
    {
        temp = x[i];
        x[i] = 0.5*x[i] + 0.5*buff[0];
        buff[0] = temp;
    }
}

EDIT2: Ok, so it turns out the problem was with the buffer size set in the library that my sound card drivers didn't like for some reason. After I changed that, everything is fine. Thank you for all the answers unfortunately I cant mark any as The Answer due to the fault being on my part. This question now can be considered as resolved.


texts:
I was reading my book where the z-transform of a signal is derived to be ${1-e^{-2bT}z^{-1}}$ . Then it goes on to say that by applying the bilinear transform we can get

$$\frac{2(1+bT+(bT-1)z^{-1})}{T(1+z^{-1})}$$

I know a little bit about bilinear transform but this example is using it in a different context. Can someone please help me understand how the bilinear transform is being used here.


texts:
I've been reading about Super Resolution Image reconstruction (Reconstruction of high resolution image from multiple low resolution aliased images contain sub pixel shifts), and i want to know why SR is possible? this document explains why SR is possible, but i didn't understand completely!

please can you explain more clearly.


texts:
If a musician (or a vocal coach) hears a human voice at a certain note, he/she can easily find the same note on piano.
Not only for human voice, if the waveform is square, or created with a synth, they will be able to figure out which note it is.
However, using DFT (or FFT) gives mixed results for any waveform other than sinusoidal. It makes sense, because we use sinusoids as base functions (or base vectors).
Now, I want to learn if there is any way we can implement a signal processing algorithm in order to detect tones regardless of their waveform. Googling didn't help. Can you help me DSP SE?


texts: My question is regarding the input size of a signal which is not a power of 2 and we have to take the fft of it. Some solutions say that suppose  if we want to take the fft of 1800 we should zero pad it till the length of 2048 to make it power of 2 and then apply the radix 2 algorithm. But there are other solutions as well which applies a combination of different algorithms without zero padding and then calculating the required FFT. My question is Does zero padding a signal to length of 2048 in case we have to take fft of 1800 makes any difference in results , if we use a combination of different algorithms to calculate the fft of size 1800. Would there be any difference or the result would be same.


texts: What is the explanation for the inverse FFT of a cross-correlated signal being shifted by some amount in the time domain?

texts:
I created two signals in the time domain: signal A which is a square wave, and signal B which is a right-angled triangle wave, as shown in the picture below.

I took the FFT of these signals, multiplied channel 1 FFT with the conjugate of channel 2 FFT. Then I inverse-FFT'd to observe it back in the time domain, and I observed this:

Now, instead if I did the correlation directly in the time domain using signal in Python like this:

txcorr = signal.correlate(wave1, wave2, mode = 'same', method='direct')

I see this output:

So, my assumption is taking the conjugate of a signal in the frequency domain reverses the time axis? Could someone give me a more mathematical background or where I could read to understand what is happening exactly in the frequency domain? Of course, when I simply do the multiplication in the frequency domain without taking the conjugate of a signal, it produces the shape expected

Additionally, this is probably the biggest confusion: Why is the "maximum point" of the cross-correlated signal in the time domain not at 0.3s? This is what I expected from what cross-correlation page in wikipedia showed:


texts: What is the true meaning of a minimum phase system? Reading the Wikipedia article and Oppenheim is some help, in that, we understand that for an LTI system, minimum phase means the inverse is causal and stable. (So that means zeros and poles are inside the unit circle), but what does "phase" and "minimum" have to do with it? Can we tell a system is minimum phase by looking at the phase response of the DFT somehow?


texts: DFT independent variable as fraction of sampling rate

texts:
When the frequency domain's independent variable is labeled as a fraction of the sampling rate, the values along the horizontal axis always run between $0$ and $0.5$, since discrete data can only contain frequencies between DC ($0$) and half the sampling rate.

But why is the highest frequency pegged at exactly $0.5$ of the sampling rate? Isn't it possible to have a scenario where the highest frequency is less than half the sampling rate?

E.g., consider a speech signal that has been filtered to remove all frequencies above $3.3$ kHz and has been sampled at $10$ kHz. The highest frequency ($3.3$ kHz) is less than $0.5$ of the sampling rate ($5$ kHz). How would the frequency domain of this signal (labeled as a fraction of the sampling rate) look like?


texts: Local signal-to-noise ratio for an image

texts:
I want to do one of these common experiments to see the effect of increased noise on the performance of a feature detector.

To do this add zero mean gaussian noise with fixed sigma to the image and look at the result. To vary the amount of noise I multiply the noise signal by different amounts before adding.

Since the detector uses a bandpass filter as a first step, to calculate the SNR of the noisy image I calculate the energy of the bandpassed original image, and the energy of the band passed noise signal (before adding them together).

Does this sound reasonable (and is that method of adding noise correct)? My concern is that because the image is being bandpassed, if I evaluate the performance of the detector at a specific point I should be using an SNR value for that local area, rather than the whole image.

How would one go about calculating local signal to noise ratio?


texts: Discrete Fourier series for a sum of deltas

texts:
so Im given a discrete sum $$x[n] = \sum\limits_{r=\infty}^{+\infty}\delta[n-rN]$$
how do I calculate its discrete Fourier series coefficients?

Thank you.

edit: this is what i've come to up to now:


texts: How are 12 outputs being reconciled so that a Homomorphic filter could be applied?

texts:
Take a look at this paper. It writes in the section 4.1:

  We selected 12 not overlapping filters, to analyze 12 different directions, rotated with respect to 15° each other. A homomorphic filter is then applied to enhance lines and to produce a dark uniform background. Finally a threshold is applied to obtain a binary mask. The threshold value is computed as mean plus standard deviation of the image intensity (see Figure 3(b)).

I could not understand this.

They began the process with a single image. Then, they applied 12 filters to that input image so that they can produce 12 outputs from one image.

How are they selecting one output image from those 12 outputs so that they can apply the homomorphic filter?  I mean, out of those 12 outputs, to which output image they would apply the homomorphic filter?

If they apply homomorphic filter and thresholding to 12 images, when did those 12 images become one image?


texts:
A question on this topic has been asked before: Combining a linear Kalman Filter with additional linear constraints.

And I checked out some of these references.

I will probably use estimate projection. However I have a couple of questions:

For the log-likelihood calculation, (i.e. when fitting the model) should I use the original or projected estimates? I presume the former ?
When forecasting should I use the original or projected estimates, in this case I think the latter ? Could someone please confirm which to use in each case and ideally why.
Lastly for the prefect measurement method shown in the link above, why is the number of restrictions forced to be less than the number of variables ? I would have thought mathematically they could also be equal to the number of variables ?

Thanks

Baz


texts: What is the difference between an equalizer and a channel estimator?

texts:
I was wondering what is the difference between channel estimator and channel equalizer?
Thanks!


texts: Probably a dumb question but I can't find the answer anywhere. It seems to me that the energy carried by a sound wave should be proportional to frequency squared. For example, the waves $x(t)=\sin(t)$, $x(t)=\sin(100\cdot t)$, and $x(t)=\frac{1}{\sqrt{2}}$ carry different amounts of energy - the last one carries no energy at all and sounds like silence. But in DSP they are all considered to have the same signal energy (integral of $x(t)$ squared). What's the significance of that?


texts: What is the meaning of the transfer function of a filter? Please explain intuitively with an example if possible.


texts: Convolution of two rectangular pulses intuition

texts:
From what I understand the convolution of the impulse response of a system with the input to that system gives the output.

Now if the impulse response is a rectangular function and the input is also a rectangular function we get a triangular function as the output.

What I dont understand is: How could this ever be the case?
As the rectangular function of the impulse response has only two levels as does the input. How can the system ever produce a value that is between these levels let alone a ramp of values between these levels.

It's obvious I am missing something but I cant quite put my finger on what that is.


texts: Detection over AWGN analysis disagrees with Monte Carlo simulation

texts:
I have a fixed vector $\mathbf{s} \in \{\pm1\}^m$ and a random vector $\mathbf{d} \in \{\pm1\}^m$.

I send either $\mathbf{x} = \mathbf{x_1} = [\mathbf{s}\space\mathbf{d}]$ or $\mathbf{x} = \mathbf{x_2} = [\mathbf{d}\space\mathbf{s}]$ over an AWGN channel to get $\mathbf{y}$.

Then $\mathbf{x}, \mathbf{y}\in \{\pm1\}^{2m}$.

From the receive vector $\mathbf{y}$, I want to detect which $\mathbf{x}$ has been sent.

Although it is suboptimal, I want to calculate the error probability of correlation detector.

The error probability of correlation detector

$$\mathbb{P}[\epsilon] = \mathbb{P}\{\textrm{declare } \mathbf{x_2} | \mathbf{x_1} \textrm{ sent}\} =\mathbb{P}\left[\sum_{i=0}^{m-1}(s_i+n_i)s_i<\sum_{i=0}^{m-1}(d_i+w_i)s_i\right] \tag{1}$$

where $s_i,d_i \in \{\pm1\}$ and $w_i, n_i \sim \mathcal{N}(0,\sigma^2)$ and they are i.i.d.

\begin{align}
\mathbb{P}[\epsilon]&=\mathbb{P}\left[\sum_{i=0}^{m-1}(s_in_i - s_iw_i)< \left(\sum_{i=0}^{m-1}d_is_i\right) -m\right]\\
&=p_\mathbf{D}(\mathbf{d}) \sum_{\mathbf{d}} \mathbb{P}_{\textrm{given } \mathbf{d}}\left[\sum_{i=0}^{m-1}(s_in_i - s_iw_i) < \left(\sum_{i=0}^{m-1}d_is_i\right) -m\right]\\
&=\frac{1}{2^m} \sum_{\mathbf{d}} \mathbb{P}_{\textrm{given } \mathbf{d}}\left[\sum_{i=0}^{m-1}(s_in_i - s_iw_i) < \left(\sum_{i=0}^{m-1}d_is_i\right) -m\right] \tag{2}
\end{align}

The random variable $\sum_{i=0}^{m-1}(s_in_i - s_iw_i)$ is Gaussian $ \mathcal{N}(0,2m\sigma^2)$. Thus

$$\mathbb{P}[\epsilon] = \frac{1}{2^m} \sum_{\mathbf{d}} \mathrm{F}_z\left( \frac{\sum_{i=0}^{m-1}d_is_i -m}{\sigma\sqrt{2m}}\right) \tag{3}$$

where $\mathrm{F}_z$ is the CDF of $z \sim \mathcal{N}(0,1)$

To check my result, I run a simulation. First I check the conditional probability "given $\mathbf{d}$"

ntest = 1000;
sigma = 1;

m = 10;
s = 2*(rand(m,1) > 0.5) - 1;
d = 2*(rand(m,1) > 0.5) - 1;

nerr = 0;
for i=1:ntest
    ys = s + sigma * randn(size(s));
    yd = d + sigma * randn(size(s));
    nerr = nerr + (yd.'*s > ys.'*s);
end

perr = nerr/ntest;  %simulated error
terr = normcdf((s.'*d - m)/sqrt(2*m)/sigma); %theory

And the values of perr and terr are similar. Good.

Next I want to check the result of (3).

ntest = 10000;
sigma = 1;

m = 10;
s = 2*(rand(m,1) > 0.5) - 1;

nerr = 0;
for i=1:ntest
    d = 2*(rand(m,1) > 0.5) - 1;
    ys = s + sigma * randn(size(s));
    yd = d + sigma * randn(size(s));
    nerr = nerr + (yd.'*s > ys.'*s);
end
perr = nerr/ntest;  %simulated error

terr = 0;
for i=0:2^(m-1)
    d = (dec2bin(i,m) - '0').';
    terr = terr + normcdf((s.'*d - m)/sqrt(2*m)/sigma);
end
terr = terr / 2^m; %theory

Theory terr = 0.0065 while simulated result perr = 0.0318. Very bad !!!

My question is that where were my mistakes ? Is my Equation (3) not correct ? Or my simulation way are wrong ? Or I have been wrong since Equation (1) and (2) ?

Thanks


texts:
I'm willing to estimate the Harmonics to Noise Ratio (HNR) of a speech signal x[k] and using autocorrelation method.

Theoretically, HNR is given as,

$$ \ HNR =  \frac{R_{xx}[T_0] }{R_{xx}[0]-R_{xx}[T_0]} $$

where $\ R_{xx}$ is the autocorrelation function (ACF) and estimated as,

$$ R_{xx,biased}[l] = \frac{1}{N} \sum_{k=l}^{N-1}{x[k]x[k-l]} $$

and,

$$ R_{xx,unbiased}[l] = \frac{1}{N-l} \sum_{k=l}^{N-1}{x[k]x[k-l]} $$

My question is how can I choose ACF function? Should it be biased or unbiased?

Constraints are signal x should be stationary (so I use the small part of the signal to be stationary) and noise is white and should be uncorreleated (N should be large enough). But these constraints has nothing to do the with type of the ACF estimation.


texts: Can multi-rate operations of decimation and interpolation can be used to implement a rate change by a factor an irrational number?

texts: I know that common rates are integer values like 2 or 3, but can multi-rate operations of decimation and interpolation can be used to implement a rate change by a factor of an irrational number? Common sense is telling me no, because irrational numbers can't be made from any set of integers by definition, so no combination of interpolation and decimation filters can generate the irrational number. There is also the fact that the MATLAB command resample.m uses the poly phase structure for implementing a sampling rate change by a rational factor. All signs are pointing to irrational numbers can't be done. Am I wrong in this thought process?


texts: FFT bin interpolation for low frequency with DC offset

texts:
I have a noisy signal in I/Q format with a peak frequency that constantly changes between 0.5Hz and 4Hz. The sample rate is 64Hz. I want an output showing the current peak frequency with a maximum lag of about 1 second.

This is my first attempt to solve this problem with FFT, and the question is: What is required to make the frequency estimate accuracy around 0.01 Hz?

The idea is to use the FFT to both do the filtering of unwanted higher frequency noise and do the estimation of the 0.5-4Hz signal peak frequency.

With one second, we have 64 complex samples go into the FFT and 64 complex magnitudes coming out. They will then represent frequencies from -32Hz to +31Hz. We look at the 6 bins representing 0Hz to 6Hz and want to interpolate here to find the peak.

The algorithm must run on a low power CPU, so computational power is quite limited and power consumption is an important parameter.

The algorithm looks like this in Scilab:

sample_rate = 64;
secs = 1;
nr = sample_rate * secs;

f_signal = 0.45; // Hz - The signals peak frequency

// Construct the ideal complex signal vector
BT = complex( cos(f_signal*2*%pi*((1:nr)-1)/nr) , sin(f_signal*2*%pi*((1:nr)-1)/nr) );

// Remove any DC offset (or not?)
// BTavg = sum(BT) / nr;
// BT = BT - BTavg;

// Make a time domain plot
figure(1);
clf();
plot(real(BT), 'red');
plot(imag(BT), 'blue');

// Perform FFT (and shift bins around to have -32Hz to +31Hz)
BF = fftshift(fft(BT));

// Create the associated frequency vector f(1) = frequency of bin 1
f = sample_rate * (-(nr)/2:(nr-2)/2)/nr;

// Make a frequency domain plot
figure(2);
clf();
plot(f(1:nr), abs(BF(1:nr)), 'blue')

// Find max peak in the first 6 bins
pmax = 0;
imax = 0;
for i=(33:38)
    p = abs(BF(i));
    printf("i = %d f = %f p = %f \n", i, f(i), p);
    if p > pmax then
        pmax = p;
        imax = i;
    end
end

// Interpolate to find better estimate of peak frequency
y1 = abs(BF(imax-1));
y2 = abs(BF(imax));
y3 = abs(BF(imax+1));
// d = (y3 - y1) / (2 * (2 * y2 - y1 - y3)); // Quadradic method
d = (y3 - y1) / (y1 + y2 + y3); // Barycentric method
fmax = f(imax) + d * (f(imax+1) - f(imax));

printf("Estimated = %.2f Hz Actual = %.2f Hz Error = %.2f Hz\n", fmax, f_signal, abs(fmax-f_signal));

The input plot looks like this:

The output from the FFT (shifted) looks like this:

The console output showing the power from the 6 bins as well as the final estimated frequency and error is here:

  i = 33 f = 0.000000 p = 44.717018
  i = 34 f = 1.000000 p = 36.588120
  i = 35 f = 2.000000 p = 12.993838
  i = 36 f = 3.000000 p = 7.911240
  i = 37 f = 4.000000 p = 5.696681
  i = 38 f = 5.000000 p = 4.459184
  Estimated = 0.24 Hz Actual = 0.45 Hz Error = 0.21 Hz

As can be seen the error is way bigger than desired for this specific test frequency of 0.45 Hz. By changing the frequency it's clear that the error near zero only for near integer frequencies.

A caveat is that the input signal also has a varying DC offset, so ideally the DC offset should be "filtered" as well.

I am sure there are ways of saving power by only calculating the 6 bins we are interested in, but that optimization only makes sense if we can get the required accuracy from the algorithm in the first place.

So again: What is required to make the error on the estimates more like 0.01 Hz?


texts: Is the DCT prone to spectral leakage like the DFT?

texts:
My understanding is that the DCT assumes the input to be an even periodic signal.
This avoids the discontinuity that can occur in the periodic extension assumed by the DFT (see image).
Does this lack of discontinuity prevent spectral leakage from occurring in the DCT?
If so, why are input signals still often windowed prior to DCT?
Thanks


texts: Dummy DSP student request for IDTFT Calculation and Integration Process

texts:
We have one DSP course in biomedical engineering, so i have been for 10 years away from mathematical calculation and working in this years by programing by c and ..., so one of our exercise is to calculating the IDTFT of some signals like this from this eBook "Digital Signal Processing Using MATLAB 3rd Edition":

  Using the definition of the inverse DTFT in (3.2), determine the
  sequences corresponding to the following DTFTs:

so the manual book solved this question by summery like this:

So, i don't know how it reached this part:

could you give me some ...?
It is possible to do this calculation by some software and seeing the formulation changing mode from base question to the answer like those shown above?

Thanks for your attention.


texts: limiter algorithm wanted: mixing two channels without clipping

texts:
I have two audio tracks (instrumental recordings) named track A, and track B. The tracks will be played in an application and the user of this application has the ability to change the balance of these tracks. The user controls a slider ranging from $x=0$ to $x=1$ with which the balance can be changed from hearing track A only, to hearing track A and B mixed equally, to hearing track B only. Mathematically, I have volume functions for track A and B, $V_A(x)$ and $V_B(x)$ respectively, that should satisfy:

$V_A(0) = 1$
$V_A(1) = 0$
$V_A(0.5)=V_B(0.5)$
$V_B(0) = 0$
$V_B(1) = 1.$

The resulting track R is calculated as follows:

$R(t) = A(t) V_A(x) + B(t) V_B(x)$ where $0 \le x \le 1.$

A naive approach would be the to let $V_A(x) = 1 - x$ and $V_B(x) = x$ but I don't want the volume of track A to go down drastically as the user moves to the equal mix at $x = 0.5$. A better approach is $V_A(x) = \sqrt{1 - x}$ and $V_B(x) = \sqrt{x}$ (see How to mix two signals without changing the overall loudness?). But ideally, I would have the following function (with the non-constant parts possibly some other power function):

$V_A(x) = \left\{
 \begin{array}{ll}
  1  & \mbox{if } x \le 0.5 \\
  2  (1 - x) & \mbox{if } x > 0
 \end{array}
\right.$

$V_B(x) = \left\{
 \begin{array}{ll}
  2x  & \mbox{if } x \le 0.5 \\
     1 & \mbox{if } x > 0.5
 \end{array}
\right.$

The problem is that both track A and track B are mastered at -0.1 dBFS (they will be played on a device and need to be as loud as possible), so obviously the last two functions will (most likely in my case) cause clipping.

I have two ideas on how to handle this:

Use my ideal function and apply some type of limiter (or would it be better to limit track A and B first, and then mix?)
Find the highest value for $V_A(0.5) = V_B(0.5)$ such that $R(t)$ never clips and then fill in the volume curve somehow. I can even limit some peaks so that this value could probably be quite high.

What are people's thoughts on this? My specific question is about idea 1. How does one go about creating a limiter? Does anyone have a good reference to some algorithms?


texts: change the frequency component of a spectrum in an audio file

texts:
I am working on an audio dynamic range compressor using Matlab.
After calculating the FFT of the audio file, I use several bandpass filters to amplify specific bands separately.
What I want to do then is remove components of a band and put them in another band.
for example, remove the components of the band 8khz-10khz and put them in the band 6khz-8khz.
would you please give me the general picture.


texts: Help me understanding the super equalizer source code

texts:
Please forgive my bad English.

Hello, I'm trying to understand Shibatch's super equalizer Winamp plugin source,
but I'm not good at DSP and mathematics.

It seems that super equalizer plugin uses overlap-add method to evaluate FIR filter.
The code that calculates the FIR filter coefficients in super equalizer is like this:
(Just the simplified pseudocode... not the exact code.)

for (int i = 0; i != FILTER_LENGTH; ++i)
{
    int t = i - FILTER_LENGTH / 2;

    result[i] =
        band_gain[0] * LPF(t, band_freq[0]) +
        band_gain[1] * (LPF(t, band_freq[1]) - LPF(a, band_freq[0])) +
        band_gain[2] * (LPF(t, band_freq[2]) - LPF(a, band_freq[1])) +
        band_gain[3] * (LPF(t, band_freq[3]) - LPF(a, band_freq[2])) +
        band_gain[4] * (IMPULSE(t) - LPF(t, band_freq[3]));

    result[i] *= WINDOW(t);
}

LPF is lowpass filter function. LPF(t, CutoffFrequency)

IMPULSE(t) is 1 when t is 0. Otherwise IMPULSE(t) is 0.

WINDOW(t) is Kaiser window function.

band_freq array contains equalizer frequency bands.

band_gain array contains the weight of each frequency bands.

I assumed there are just four frequency bands.

band_gain[0] = weight of 0Hz ~ band_freq[0]Hz
band_gain[1] = weight of band_freq[0]Hz ~ band_freq[1]Hz
band_gain[2] = weight of band_freq[1]Hz ~ band_freq[2]Hz
band_gain[3] = weight of band_freq[2]Hz ~ band_freq[3]Hz
band_gain[4] = weight of band_freq[3]Hz ~ SAMPLING_RATE / 2

I want to know this part:

result[i] =
    band_gain[0] * LPF(t, band_freq[0]) +
    band_gain[1] * (LPF(t, band_freq[1]) - LPF(a, band_freq[0])) +
    band_gain[2] * (LPF(t, band_freq[2]) - LPF(a, band_freq[1])) +
    band_gain[3] * (LPF(t, band_freq[3]) - LPF(a, band_freq[2])) +
    band_gain[4] * (IMPULSE(t) - LPF(t, band_freq[3]));

What formula does the code using?
Or... what should I study to understand the code?

Any answer is welcome.
Keyword suggestion is welcome.

Thanks.


texts: How can one improve the robustness of adaptive beamformers to signal mismatches?

texts:
Compared to standard (conventional) beamformers, adaptive beamformers such as minimum variance distortionless response (MVDR) have much better resolution and interference nulling capabilities. Briefly, the MVDR beamformer is as follows:

If $\mathbf{a}_\theta$ is the look-vector in the direction $\theta$ and if $\widehat{\mathbf{R}}$ is the covariance matrix of observations, then the adaptive weight vectors are calculated as

$$\mathbf{w}_\theta=\frac{\widehat{\mathbf{R}}^{-1} \mathbf{a}_\theta}{\mathbf{a}^H_\theta \widehat{\mathbf{R}}^{-1} \mathbf{a}_\theta}$$

and the output power of the MVDR beamformer is

$$B_{mvdr}(\theta)=\mathbf{w}_\theta^H \widehat{\mathbf{R}}\mathbf{w}_\theta$$

For the sake of completeness, the conventional beamformer output power is $B_{conv}(\theta)=\mathbf{a}_\theta^H \widehat{\mathbf{R}}\mathbf{a}_\theta$.

However, these adaptive beamformers are also prone to severe performance hits when there is a signal mismatch i.e., when the steering vector model for the signal of interest doesn't match what's in the data or when the actual array positions are perturbed from the array positions used in the steering/look vectors (which is what happens in practice). The performance degradation can be so large that it becomes inferior to the conventional beamformer.

How can one increase the robustness of adaptive beamformers (in this case, the MVDR, but applicable more generally) to signal mismatches?

Example:

Here's the result from a simulation with two plane wave signals, one at $+30^\circ$ and the other at $-22.5^\circ$ and 10 dB down from the first, impinging on a 10 element uniform line array (spaced at $\lambda/2$, where $\lambda$ is the wavelength) in a noisy environment (the exact value of the noise power is immaterial here and you can assume it to be much lower than the signal powers). The covariance matrix is formed from 50 observations (snapshots). The top two panels show the outputs of the MVDR and the conventional beamformers and you can clearly see that the MVDR has high resolution compared to the conventional. Both do a good job of estimating the signal power.

However, when the element positions are slightly perturbed (perturbation is unknown to the beamformer), causing a mismatch in the actual signal steering vector and the model, the signal power estimates of the MVDR beamformer goes kaput. It can hardly measure them correctly, although it still does a good job of localising the signals. Comparatively, the conventional beamformer gets the relative powers right, although its harder to localise the signal due to increased sidelobes. (Note: The plots have all been normalized by the maximum power so although you can't see it, the peak at $30^\circ$ in the top left plot is at 0 dB. So this doesn't truly show how the actual power estimate degrades due to mismatches, as I just wanted to paint a qualitative picture.)

http://dsp.stackexchange.com/questions/146/how-can-one-improve-the-robustness-of-adaptive-beamformers-to-signal-mismatches#comment-308 :

  The title for the lower left plot should read MVDR: Perturbed array.
  I'll correct the mistake sometime soon.
  – yoda Sep 3 at 5:28
  -->

Is it possible for me to improve the estimates of the signal powers using the MVDR beamformer in the perturbed case?


texts: monochrome (1-bit black-and-white) image conversion

texts:
When converting scanned text to a 1-bit black-and-white image what are some filters that one can apply in the process to improve the result? Right now I am running into the problem where dithering errors make the image look horrible.

Update:
I think trying to undo dithering is a much harder problem. How can I convert the first image to a monochrome image? The default approach convert -monochrome img1 img2 is shown below. I've also tried a two-step approach: 1) decrease depth (color palette) and gamma 2) convert to a bilevel image (not shown). Other things I played with included imagemagick's ordered-dither (at various settings), but it wasn't as good as the two-step approach.


texts: Phase response distortion due to IQ demodulation

texts:
I have a bandpass signal, I do IQ demodulation meaning:

first step:

$IQ(t) = x(t)*exp(-2\pi f_ct)$

second step:

low pass using $filtfilt$ function in MATLAB. (in order to avoid phase distortion)

Here, $F_c=5 MHz$ and $F_S=40 MHz$

For low-pass filter I use this:

[b,a] = butter(5,2*Fc/Fs) ;
IQ(t) = filtfilt(b,a,IQ(t))*2;

The plot shows frequency response of $x(t)$, $real(IQ(t))$ a.k.a $I$, $imag(IQ(t))$ a.k.a $Q$ and $IQ(t)$, respectively.

When I reconstruct this signal again to the original bandpass signal:

$x'(t) = IQ(t)*exp(+2\pi f_ct)$

The magnitude, phase response and power spectrum of the original and reconstructed signal is a follows:

note: To plot phase response I use matlab's $unwrap$ function.

The frequency spectrum of original and reconstructed signal:

The phase response of the $x'(t)$ seems distorted at $2f_c = 10MHz$ .

Question1: Is this normal?

Question2: Does this happen because low-pass filter is not reversible?

Question3: Can I recover that phase response?


texts: Frequency response and sampling theorem for triangular function

texts:
The triangular function is defined as follows:
$h_l(x) = \begin{cases}1-|x|,&|x|<1;\\0&\text{otherwise}.\end{cases}$
According to ccrma.stanford.edu:
"If the output of the interpolator is also sampled, this can be modeled by sampling the continuous-time interpolation result, thereby aliasing the $\text{sinc}^2$ frequency response. [...] The Fourier transform of $h_l(nT/L)$ is the same function aliased on a block of size $ f_s=L/T$ Hz."
What I am thinking: If the Fourier transform $\hat{f}$ has bounded support (bandlimited), then $f$ has unbounded support (infinite length). This is a consequence of the uncertainty principle. But the triangular function (i.e. $f$) is $0$ for $|x| > 1$, so $f$ has bounded support (and $\hat{f}$ cannot be bandlimited). The sampling theorem requires $\hat{f}$ to be supported on some interval $[-B, B]$ (bandlimited). Then we cannot apply the sampling theorem, because $\hat{f}$ is not bandlimited. So aliasing will occur because we cannot find a good sampling rate.
Is my interpretation correct? Aliasing will always occur for $h_l$ no matter whether we look at $\text{DTFT}(\text{sample}(h_l))$ or at $\hat{h_l}$ (Fourier transform)? Then why is linear interpolation so popular in image processing (bilinear interpolation) if it has such a bad frequency response?


texts: Low analog lowpass + ADC vs High analog lowpass + ADC + low digital lowpass

texts: I'm doing EEG data collection and my amp forces use of an analog lowpass filter before the ADC. I'm fairly sure that I won't be interested in phenomena associated with frequencies of 30Hz or higher, so it strikes me that the most disk-space efficient thing to do is to analog lowpass at 30Hz (it uses a Bessel filter) and tell the ADC to sample at 60Hz. However, disk space considerations aside, I'm curious whether there is any difference between this strategy and, say, setting the lowpass to something really high like 250Hz and sampling at 500Hz, then once digitized and doing analysis, applying a digital lowpass at 30Hz. Would the double filtering induce artifacts? (i.e. Should I opt for the low analog lowpass for reasons other than disk space?)


texts: Using a Kalman filter to predict samples between measurements

texts:
(cross posted from stack overflow)
I have a dataset where measurements were taken at 10 Hz, and I am trying to use a Kalman filter to add predicted samples in between the measurements, so that my output is at 100 Hz. I have it working ok when the velocity is linear, but when the direction changes, the filter takes a while to catch up. I am new to Kalman models, so am very likely making some mistakes in my settings. See image

for an example, the red is measured data, with stepping in between measurements. The blue is the Kalman corrected.

std::vector measurements is a dummy data array I am testing with.

The main Kalman code is based on this:

https://github.com/hmartiro/kalman-cpp/blob/master/kalman.cpp

Am I correct in feeding the prediction back into the filter for every loop? Or is this line:

yPos << kf.state().transpose(); wrong?

my code is:

int main(int argc, char* argv[]) {

  int n = 3; // Number of states
  int m = 1; // Number of measurements

  double dt = 1.0/30; // Time step

  Eigen::MatrixXd matA(n, n); // System dynamics matrix
  Eigen::MatrixXd matC(m, n); // Output matrix
  Eigen::MatrixXd matQ(n, n); // Process noise covariance
  Eigen::MatrixXd matR(m, m); // Measurement noise covariance
  Eigen::MatrixXd matP(n, n); // Estimate error covariance

  // Discrete motion, measuring position only
  matA << 1, dt, 0, 0, 1, dt, 0, 0, 1;
  matC << 1, 0, 0;

  // Reasonable covariance matrices
  matQ << 0.001, 0.001, .0, 0.001, 0.001, .0, .0, .0, .0;
  matR << 0.03;
  matP << .1, .1, .1, .1, 10000, 10, .1, 10, 100;

  // Construct the filter
  KalmanFilter kf(dt,matA, matC, matQ, matR, matP);

  // List of noisy position measurements (yPos)
  std::vector<double> measurements = {
     10,11,13,13.5,14,15.2,15.6,16,18,22,20,21,19,18,17,16,17.5,19,21,22,23,25,26,25,24,21,20,18,16
  };

  // Best guess of initial states
  Eigen::VectorXd x0(n);
  x0 << measurements[0], 0, 0;
  kf.init(dt,x0);

  // Feed measurements into filter, output estimated states
  double t = 0;
  Eigen::VectorXd y(m);

  for(int i = 0; i < measurements.size(); i++) { //ACTUAL MEASURED SAMPLE

      yPos << measurements[i];

      kf.update(yPos);

      for (int ji = 0; ji < 10; ji++)  // TEN PREDICTED SAMPLES
      {
          t += dt;

         kf.update(yPos);

          yPos << kf.state().transpose(); //USE PREDICTION AS NEW SAMPLE

      }
  }

  return 0;
}

Thank you.


texts:
I am a total newbie in signal processing.
I have been trying to plot on simulink a time/frequency domain noise in PLC (Powerline Communication) that matches a real-world research "Time Frequency Analysis of Noise Generated by Electrical Loads in PLC", however, I failed to plot such noise.

Generally speaking about my aim, I am trying to collect appliances’ noise that causes bit errors in PLC by comparing the original bit sequence with the received bit sequence. So, I need to simulate each appliance noise individually like the article mentioned above. Unfortunately, I couldn’t find any data source where I could download appliances’ noise, so I am creating my own simulation to collect my appliances’ bit-errors data.
I have been able to simulate PLC using FSK modulation/demodulation with 1000 bits per second as a binary sample rate for my bit-sequence. For the FSK modulation, I have 7500 Hz as a carrier frequency for the bit-sequence to the AC channel that has 50 Hz/220 V as the typical AC source voltage, but now I am stuck on the phase of creating Appliances’ noise which has to matches real-world noise.

If anyone knows how to plot on Simulink as I don't know how to plot such signals:

In frequency domain: background noise as high as 55dB for the measurement frequency range with the power contained up to 4MHz
In time domain: Raise up to 45dB above the background noise floor for most of the AC cycle duration

As far as the article shows, it should look like this.

I need to plot a huge dataset, so if anyone could help me understand how to plot such noise, I will be able to plot the rest.


texts: audio processing for audio files quick listening

texts:
I need to speed the listening task. I want to listen only 1 second every three seconds.
Could you give me some advice?
Many thanks!
Best regards, Sergio


texts: How to best separate my signal into frequency regions?

texts:
I have a signal with frequencies betwwen [0.5,500] Hz and i want to create a filter to best separate the signal in the following regions:

[0.5-3]Hz
[3-8]Hz
[8-13]Hz
[13-30]Hz
[30-500]Hz

My questions are:

How can i find the sampling frequency if i know that the signals are
'downloaded' in Europe (energy transport network of 50Hz)
Whats the most appropriate computational method to separate the signal in
the regions.

I think that in order to find the sampling frequency i need to find the max frequency of my signal and then fsampling=2*fmax.Does the information about the 50Hz network affects my sampling frequency?


texts: Mahalanobis Distance For Pairwise Distance

texts:
First post here so please be easy on me.

I have an application where I have one data point in high dimensions, represented by an $m \times n$ matrix. Call it $X_1$. And I have to compare this data point to many other data points of the same dimension, call them $X_k$ where $k\in \{2,3...,z\}$, give each comparison a score, and then choose the $X_k$ that is most similar ( closest to ) to $X_1$. The way I am currently implementing this score is just the Euclidean distance (Frobenius Norm). Let $s_k$ be the score,

$$s_k = ||X_1-X_k||_F$$

So now I am taking a class where we are looking at many different distance measures. I know my data set has some covariance structure in it, not IID, so it has come up that perhaps a better way of looking at distance is with the Mahalanobis distance, which normalizes the distance from the mean using the covariance matrix, so that distances are relative their positions along the principal components. Here $\hat{x}$ and $\hat{\mu}$ are vectors representing the data point and the mean of the distribution.

$$D_M=\sqrt{(\hat{x}-\hat{\mu})C(\hat{x}-\hat{\mu})^T}$$

My question here is: can I use the theory of the Mahalanobis distance with my scoring function instead of using the Frobenius norm?

The issues I see are that the Mahalanobis distance is not measuring a pairwise but instead the distance from a point to the mean of the distribution. Another issue I see is that my data points are 2-dimensional matrices and the Mahalanobis distance is implemented on vectors.

Thanks


texts: Bandpass filter on noisy sinusoid that predicts turning points

texts:
Is there an application where it is required to predict turning-points of a signal that
is like a sinusoid with constant frequency, variable amplitude and some noise (Gaussian additive noise)? The amplitude of the noise should be less than 10% of the amplitude of the original sinusoid every "oscillation".

Reason I ask is that I think I found an interesting causal filter for this very specific
situation, but I have no imagination for the moment on how it could be usefull.

My causal band-pass filter turns earlier than the turning point of the clean sinusoid,
so it's like the filter has a negative lag. And the computation time is short, less than my code using FFT (in R: statistics software).

I'm a hobiest in DSP so please excuse my ignorance. My knowledge about filters comes mainly
from finance (econometrics).

To the man with a hammer, the entire world looks like a nail..

Thanks for any suggestions.


texts:
How to prove a cosine signal with irrational period and another cosine signal with a rational period are orthogonal?

For example, $\cos(4t)$ and $\cos(4\pi t)$ sum is aperiodic. Hence it is said that these two signals are orthogonal over infinite limits? But how to prove this? Normal integration of product of the two signals cannot be used because we get $\cos(\infty)$ which is unknown.


texts: How to distinguish between recorded and a live vocal signals?

texts:
Is it possible to detect if a sample of human voice by smartphone/computer(44Khz) is sourced from another electrical device or originated  from a human being

I want to build a model in matlab such that, samples a vocal signal and decides whether it is a real human being talks a live or recorded, is it even possible?


texts: Digital implementation of type I compensation network

texts:
I'm looking for advice or references on how to implement a "type I" compensation network of digitally.  An analog example is shown in the following diagram:

I am familiar with the basics of digital signal processing, and have read about digital integrators.  Would it be possible to use something like a first order IIR digital integrator, and just flip the sign of the output depending upon whether the input is above or below VREF? The goal is to drive a PWM output of a microprocessor, so there may be a better way to do this in the digital domain than attempting to copy the analog implementation.  Any advice would be appreciated.


texts: Best averaging to convert 50Hz Chi data to one hour averages?

texts:
I have a $50\text{Hz}$ data set of temperature variance destruction or $\chi$ in the ocean for 365 days and hence is too bulky to handle. I want to convert it to hourly averages.

Simple question: How do I decide the best averaging method given that I do not have any pre-held notion about the data i.e I do not know what to expect? How can I justify its use for the current data? How do I make sure I do not lose the seasonal, diurnal trends etc. when I average?

My further calculations heavily depend on this processing and also this processed signal will give me spectra of different physics in the ocean. So I want to tread carefully. Any literature, suggestion should do to start off.

PS: In few posts and forums, I have read people commenting on moving averages as a "lousy" way to do it.


texts: Changing the number of cycles of an sinusoidal input signal

texts:
Imagine we have an ultrasonic sensor that transmits a multiple cycles of sinusoidal signal/waveform that can travel up to a specific distance/range.

Can we achieve more distance coverage using more number of cycles in the input signal (at the same frequency) or this change can not affect the coverage?

I understand that more number of cycles will result in more signal energy. Does that mean that having more energy in the signal will result in less signal deterioration when SNR goes down in longer distances due to attenuation?


texts: From fft result one can affirm that the signal is affected by white noise?

texts:
I calculated the FFT of the Az acceleration sampled with MPU6050 connected to Arduino.

Wikipedia definition: "White noise" is  a random signal having equal intensity at different frequencies, giving it a constant power spectral density

Can I say from my FFT that the reason of the variance of the signal, in my test condition, is caused by white noise?

Additional question on FFT

If my sampling frequency is 56Hz are my bin_frequency 28?
My frequency bin is influenced by the length of the input vector or only by my sampling frequency?

Test condition:

I recorded of the z value in static condition on a table, without any perturbation
The signal was  nonuniformly sampled so I resampled it

Input

Sampling frequency 55.556

Az 718X1 Vector

Std_Input_Vector_Resampled 0,0385

FFT Output

FFT 718X1 Vector

Without Zoom

ZoomX3

ZoomX4

This is my first time I approach the FFT(I've been studying this argument from one month so I am very noob), and I have a civil engineer background, so is my first time with signal analysis, never studied Signal Theory before.

My interest for the FFT is to define the best filter to identify the pothole and road anomalies

Thanks for your time, thanks for your patience.


texts: Help with a product of a sum of infinite series

texts:
i would like some help with this:

given:
$$
h[n]=2^{-n} u[n]
$$
calculate the series:
$$
 \sum_{r=- \infty }^ \infty  h[n+4r]
$$


texts: What is the transfer function of this system y(n) =2y(n-1)? Here there is no input dependence.  If you give one output sample, all the future output samples are determined by multiplying with 2.


texts: Why the spectral coherence is unity for all frequencies between single-frequency time series and itself

texts:
In the example below, I am plotting the coherence between time series and itself. The time series do has one frequency.The coherence magnitude was one for all frequencies. I wonder why it is not zero at all frequencies except at the time series frequency. I found that if I include random variable, I get zero coherence at most frequencies except the time series and its surrounding frequencies.
Here is my questions:

Why the coherence is unit for all frequencies however the time series have only single frequency.
Why I got more realistic (coherence values at the time series frequency, and zero else where) values when I add noise to the time series.
Why the coherence (with the noise added case) show more unity for more that one frequency. It looks like that the coherence have unit at the time series frequency and the surrounding frequencies.

import matplotlib.pyplot as plt
from scipy import signal

nt=10000
w1=0.02

xp_pre=np.zeros((nt))
yp_pre=np.zeros((nt))

for t in np.arange(nt):
    xp_pre[t]=np.cos(2*np.pi*w1*t)
    yp_pre[t]=np.cos(2*np.pi*w1*t)

#xp_pre=xp_pre+np.random.randn(nt)*0.03
#yp_pre=yp_pre+np.random.randn(nt)*0.03

plt.figure();plt.plot(xp_pre);plt.plot(yp_pre);plt.show()

#f,coh = signal.coherence(xp_pre,yp_pre,window=signal.get_window('hanning',200),noverlap=100)
f,coh = signal.coherence(xp_pre,yp_pre)
plt.figure();plt.plot(f,coh,'*')
plt.show()

Thanks


texts: In-Phase and Quadrature components of low-pass equivalent of a band-pass signal?

texts:
I was wondering what is the physical meaning of IN-PHASE and QUADRATURE components?
In the context of low-pass equivalent of band-pass signals.
Any help is appreciated!


texts: What does the Phase part of a FFT of a received signal represent?

texts:
Imagine we've receive a signal, $R$, at a distance of $L$ from a transmitter.
Then, we do FFT on the received signal. This is what FFT gives us:

\begin{align}
R(f) &= A e^{j\phi}   \tag1\\
\
\end{align}

where, $A$ is the amplitude, we can take it by abs() function over the FFT. And, $\phi$ is the phase and we can take it by angle() function over FFT.

Now this is my confusion:
If the initial phase, I mean the phase of the transmitted signal, would be $\phi_0$, then,
Is it correct if I say:

$\phi = \phi_0 + \theta$ ?

$\theta$  is a number that is added to the transmitted phase over the link.

or I can also write this:

$\phi_{received} = \phi + \phi_0$ ?

$\phi_{received}$ is the phase of the received signal.

Which one of those is correct?

Please take this into account that I have no information about the transmitted signal. Just assume a transmitted signal with initial phase $\phi_0$.

Also just consider a very simple case, like a line-of-sight signal.
And, I'm using MATLAB for FFT, and angle().

EDIT:
I have two received signals at the distances of $L_1$ and $L_2$, with the phases of $\phi_1$ and $\phi_2$.
I calculated the difference between these two phases. And, plotted them.

And, the phase shift between these two received signal should be

$\phi_2 - \phi_1 = 2\pi f \tau$

So, Can I just go and consider the two signals in the time domain and calculate the time difference, $\tau$, between them?
Any thoughts?
What am I missing here?


texts: Phase information from product of complex conjugate transfer functions

texts:
For an input signal $x(t)$ and output signal $y(t)$ through an LTI system I $H(t)$ I have found and interesting property that concerns signals' power spectral densities:

$${\lvert H(\omega)\rvert}^2 = \frac{\phi_y(\omega)}{\phi_x(\omega)}$$

This relation gives me information on the magnitude of $H(\omega)$ but no phase information. Is there a way to obtain phase information of $H(\omega)$ from this relation?


texts:
I'm reading about fault detection via signal processing in time domain. One possibility is to check that first derivative of the signal is in some predefined bounds. The text says that to obtain the first derivative of the output signal $y(t)$, I can use this analog filter:
$$F(p)=\frac{p}{\frac{1}{\tau}p+1}$$
for some sufficiently large $\tau$. Why is that? I don't  know much about analog filters. I only know that Laplace transform of the first derivative of a function $f(t)$ with zero initial conditions is $pF(p)$. I also understand that $X(p)=\frac{1}{\frac{1}{\tau}p+1}$ is Laplace transform of exponential decay function $e^{-\tau t}$. Could you please explain to me why the filter needs to be designed like that and how does it work physically?

Thank you.

My wrong reasoning: If $u(t)$ is the input signal and $y(t)$ is the rate of that signal, then
$$
y(t) = u(t)'. \\
$$
The transfer function would then be
$$
\begin{align}
Y &= pU \\
F &= \frac{Y}{U} = p.
\end{align}
$$


texts:
I'm trying to use a cardinal sine as a lowpass filter for a cosine signal with a fundamental of 1k.

In frequency domain, what really happens is that I'm multiplying two impulses (centered at -1k and 1k) with a rectangular pulse of width equals 2 (convolution is multiplication in frecuency domain). The result would be a constant zero function. However using the above code in Matlab I'm getting again the sinc function as output. Any help would be appreciated (attached the full plot).

D = 5;
f0 = 1000;
fm = 2*f0;    % nyquist reestriction

t = (-D:1/fm:D);
x = cos(1000*2*pi*t);
h = sinc(t);

x_p = filter(h,1,x);
plot(x_p);


texts:
At first, I thought harmonics come from the signal being periodic. However, we know a sine wave is also periodic but contains a single frequency and no harmonics.
A harmonic for me is a component (pure sine-wave) with a frequency $f_k = k f_0$ for some $k\in \mathbb{Z}$ and a certain frequency $f_0$ called the fundamental frequency. To this, I add the condition that its amplitude is not zero.
We can also imagine the sum of two different frequencies, one is not necessarily the harmonic of the other, and get a periodic signal but no harmonics (according to my definition, which might need correction if it's wrong), just the two frequencies.
Below, the result of a simulation where I have the sum of two sine-waves with frequencies 100 Hz and 201 Hz, respectively:
$$
x(t) = \sin(2 \pi 100 t) + \sin(2 \pi 201 t)
$$
The signal is periodic (with a period equal to 1 s) but it does not contain harmonics (again according to my definition!).

Now, I'm wondering what's the real root-cause for harmonics appearing on a signal.


texts: Spectral baseline distortion simulation - Bruker smile

texts:
I work with Magnetic Resonance Spectroscopy data. Data (complex signals) obtained from Bruker scanners (one of the biggest manufacturers of MRS scanners) start with an important group delay - around 60-80 points (in the time domain) - probably due to the fact that they use a long FIR filter.

The effect of this group delay cumulated with some other (unfortunately unknown) problems results in a baseline distortion called Bruker smile/frown. (in the frequency domain, after FFT)

I would like to simulate this kind of distortion (smile or frown).
I tried with a bandpass FIR of several taps (80- 140) and I cannot reach this effect.
Could you recommend me any way of getting this distortion? Thank you in advance.


texts: Reconstructing band-limited signal's integral from discrete samples

texts:
I'm thinking about building an e-bike consumed power meter.

I'm planning on sampling a low-pass filtered current sensing and voltage analog signals on the battery to e-bike controller connection, then somehow digitally integrating the product of these 2 (instantaneous power) to get the consumed power.

I could over-sample say by 10 or 20X of the analog bandwidth and then use a simple sum of areas of the trapezoids to get a reasonably accurate integral (consumed power).

But I'm thinking maybe there's a relatively simple DSP technique to get the integral without such significant oversampling that can be implemented in an FPGA? (e.g. to get a 5% accuracy with sampling at 3x of the analog bandwidth)

Some digital low-pass filter and then the same sum of trapezoid areas on the result?


texts:
I am almost done with design and implementation of my first FIR filter using C++. Its a Sinc-Windowed FIR filter designed for a cut off frequency of 100 Hz with sampling rate of 500, with a filter order M = 100. All the operations have been carried out in discrete time domain. So basically my question is regarding testing to see if my filter is generating the desired response .

I presume one way to test the output would be to apply DFT to it and check the response spectrum of the filter in the frequency domain, but for this I have to develop a function that performs the DFT on the output .

Question :To avoid this, I am looking,

---->Are there any methods that make it possible to test the filter response in the (without the need for applying DFT) time domain itself?.

----> If the above is not possible, are there are any alternative methods like some C++ DSP library with FT function's you can suggest?

Or Any suggestions and advice regarding this with your experience would be appreciated. Thanks in advance!


texts: Is it a correct way to calculate the Mean Square Displacement as function of time?

texts:
I would like to know if I calculate the MSD correctly, here is a Jupyter snippet I am using:

import numpy as np
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
%matplotlib inline

x = np.random.randn(100)

df = pd.DataFrame(x, columns=['x'])
df['time'] = np.arange(0,len(df))
df['MSD'] = pd.Series((x - x[0])**2).rolling(window=10, min_periods=1).mean()

x_plt, = plt.plot(x, label='x(t)')
msd_plt, = plt.plot(df['MSD'], label='MSD(t)')
plt.legend(handles=[line_up, line_down])

Is there a smart way to average over the variances? The original definition of MSD has average inside, but I don't relay averaging here because I have only one particle.


texts: In Search of Voice Activity Detection Algorithms - What happens to the frequency or pitch or speech features as speaker moves away from microphone

texts:
Actually i am in  search of a Voice Activity Detection Algorithm which could distinguish between voice and non-voice

Roughly speaking it must not detect even a bullet sound,even a foot stepping and other non speech activity should only detect people conversation  or any one shouting

in that search this question arises in my mind and i want know the effect of noise and distance between speaker and microphone on the speech features like pitch,frequency,cepstrum,zerocrossing rate,power spectral density,entropy  etc

if some component wont get distorted i would like to extract that feature and do the activity decision on that

Can any one help me in extracting dominant parameter of speech which would differentiate it from other common sounds even in Lower SNR conditions <0dB

Note:my algorithm  expects voice activity happens at a distance of at least 10m away from microphone and continuous generator hum as background noise


texts:
Why exactly is  X(0) the DC component of a signal?

How is it equal to N times x(n)'s average value and why it is at X(0)?


texts:
$$\hat f(x, y) = g(x, y)-\frac{\sigma_n^2}{\sigma_L^2}\left[g(x, y)-m_L\right]$$

What are the meanings of the following terms:

$m_L$
${\sigma_\eta}^2$
${\sigma_L}^2$

Here we see that $m_L$ is subtracted from the Image and then the whole term is multiplied by $\frac{{\sigma_\eta}^2}{{\sigma_L}^2} $. Then, the whole term is again subtracted from $g(x,y)$.

What does that actually mean?


texts: Multisampling of an ADC channel in a PWM cycle in TI C2000

texts:
I want to sample an ADC channel multiple times in one PWM cycle.
I can trig any ADC channel in anywhere in my PWM timer counter now. But I want to trig and sample same channel in many times within a pwm period.
I tried to explain by drawing below.

Is there any solution for this?
Thanks and best regards


texts: How to compute the energy of two signals?

texts:
Say I have two signals, one being $g_1(t)$ and the other one being $g_2(t)$, and say that I am being asked to compute the energy of $g_3(t)$, where $$g_3(t) = g_1(t) + g_2(t)$$

Is it correct to compute the energies of $g_1(t)$ and $g_2(t)$ separately and simply add the two values to find $g_3(t)$ or is it correct to compute them as one value?  I.e. $$\int_ {-T}^{T} \big[g_1(t)+g_2(t)\big]^2 \,dt $$


texts: Spectrogram of a single tone complex signal has two dark lines?

texts:
I am trying to plot the spectrogram of a complex signal I generated. I have written code to generate this signal and plot the spectrogram. It works. However I see that there are two dark lines present within the spectrogram. I'm curious to know the reason for the existence of these. I have also noticed that by varying the NFFT, the dark lines tend wobble around.

I paste my code below:

import matplotlib.pyplot as plt
import numpy as np
import scipy.signal as signal

def plot_spectrogram(data,NFFT,Fs,ex ):
    plt.specgram(data, NFFT=NFFT, Fs=Fs)
    plt.title("Spectrogram of data")
    plt.ylim(-Fs/2, Fs/2)
    plt.show()
    plt.close()

    if ex:
        exit()

###  Parameters
F_offset = 250000
Fs = 1140000

###  Generate a digital complex exponential with phase -F_offset/Fs
fc1 = np.exp(-1.0j*2.0*np.pi* F_offset/Fs*np.arange(len(x1)) )
plot_spectrogram(fc1, 512, Fs,ex=True)


texts: Problem with scipy.signal.butter for designing filter bank

texts:
I need to design a filter bank of band pass filters (code for implementation is shown below). There is a problem in that if my center frequencies fc are much smaller than the sampling frequency fs then you get the behavior as shown in the first graph:

A temporary solution is to change the sampling frequency to fs = fh*4 where 'fh' is the hi cutoff frequency of each respective band. This gives

Which is what I want, however because the sampling frequency is no longer the true one (relative to the time series I then want to filter) I cant use it for filtering purposes using lfilter(b,a,x) where x is my time series because all of the bands will be identical since they have the same non-dimensional frequency! My question is what can I do to make this work as expected for a true fs? Presumably there is an issue with finite precision?

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import butter, lfilter
from scipy.signal import freqz

def butter_bandpass(lowcut, highcut, fs, order=9):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a

def fbank(fs):
    """
    Calculate third octave filter bank.
    An octave is the interval between two frequencies having
    a ratio of 2:1. Hence for 1/3 octave
    fh = fc*(sqrt(2)^(1/3))
    fl = fc/(sqrt(2)^(1/3))
    """

    # ISO standard center frequencies (bins)
    F0 = [50, 63, 80, 100, 125, 160, 200, 250, 315, 400, 500,
          630, 800, 1000, 1250, 1600, 2000, 2500, 3150, 4000,
          5000, 6300, 8000, 10000, 12500, 16000, 20000]

    plt.figure(figsize=(10,6))
    for fc in F0:
        fh = fc*(2**(1/6.0))
        fl = fc/(2**(1/6.0))
        fs = fh*4 # Temporary solution
        fs = 5e4
        b, a = butter_bandpass(fl,fh,fs)
        w, h = freqz(b, a, worN=2000)
        plt.loglog((fs * 0.5 / np.pi) * w, abs(h), label="fc = %d" % fc)

    plt.xlabel("Frequency (Hz)")
    plt.ylabel("Gain")
    plt.legend(loc='best')
    plt.savefig("test_filter.png", bbox_inches='tight')
    plt.show()
    return 0


texts:
Could you help me with this theoretical question.

The purpose is to design a discrete-time system. What it does is, it finds the maximum value in a one-dimensional sequence whose elements consist only of "non-negative" real numbers and contains only a single global maximum with no local maxima.

Thanks in advance.


texts: Is the discrete Gaussian kernel an eigenfunction of the DFT?

texts:
So the Gaussian function is an eigenfunction of the Fourier transform because it transforms into itself, right?
But this isn't true for the sampled Gaussian in the DFT because the tails of the function are truncated, right?
Wikipedia describes a discrete Gaussian kernel here and here (solid lines), which is different from the discretely-sampled Gaussian (dashed lines):

the discrete counterpart of the continuous Gaussian in that it is the solution to the discrete diffusion equation (discrete space, continuous time), just as the continuous Gaussian is the solution to the continuous diffusion equation

Does that mean it also DFT transforms exactly into itself?  If not, is there a similar Gaussian-like function that does?


texts: Frequency of a noisy signal with fixed width pulses

texts:
My electricity meter has a front panel LED which pulses in proportional to the current energy rate. The frequency changes but the pulse width is set at 250ms.  Due to environmental conditions there's plenty of noise. Currently I'm sampling this at 10kHz, the pulses may be a few seconds apart.

I would like to filter these pulses and ignore anything else that wasn't 250ms long (ie the noise) and find the frequency of the train. Are there any suitable methods to achieve this?

My first thought was to differentiate the signal but this would be easily confused by large noise spikes.
I also wondered about testing it's coherence with a simulated pulse train, but as the frequency changes this would be potentially infinitely intensive AFAIK.

If I was to build an electronic solution I'd use edge comparators to find the pulses.  I'm wanting to use NumPy/matplotlib and deal with the signal in the frequency domain, if possible (purely as an educational exercise).


texts: What is the correct solution for Fourier transform of unit step signal?

texts:
The unit step signal defined as

$$
 u[n]= \lbrace 1; n>=0; \\
               \qquad0; n<0 \rbrace
$$

has three possible solutions for its Fourier domain representation depending on the type of approach. These are as follows -

The widely followed approach (Oppenheim Textbook)- calculating the Fourier transform of the unit step function from the Fourier transform of the signum function.

$$
     F(u[n])=U(j\omega) = \pi\delta(\omega)+\frac{1}{j\omega}
$$

Fourier Transform calculated from the Z transform of the unit step function (Refer Proakis Textbook, Digital Signal Processing Algorithms and applications,  pages 267,268 section 4.2.8 )

$$
     U(j\omega) = \frac {e^{\frac{j\omega}{2}}}{2j\sin \frac{\omega}{2}};   \omega \neq 2\pi k; k=0,1,2,3...
$$

Fourier transform calculated by splitting into even and odd functions - followed in Proakis Textbook (Refer Proakis Textbook, Digital Signal Processing Algorithms and applications, page 618 section 8.1 )
$$
      U(j\omega) = \pi\delta(\omega)+\frac{1}{1-e^{-j\omega}}
$$

The 2nd representation can be ignored since it is not a well-behaved function. But the approaches followed by Proakis and Oppenheim are equally valid (they extend the Fourier transform to include impulses in the frequency domain)
But the confusion is that they provide different solutions.

Is there any mistake in my understanding? or am I missing any crucial point?? Kindly help me understand this and the correct form that can be used in all applications. (I found that the Oppenheim approach is used in deriving the Kramers-Kronig Relations and the Proakis approach used in the derivation of the Hilbert transform)


texts: Audio watermark embedding text attack resistant

texts:
I'm developing an audio watermarking algorithm which embeds text (binary form) in DWT-DCT coefficients. The goal of the algorithm is to hide song information (like artist, song, album,...). Everything works fine the text is recovered using another algorithm which extracts the hidden text. However my algorithm doesn't resist attacks, I use Stirmark for Audio v02 which I found here. The quality of the audio is not really affected, in fact I get good values when using the Normalized Cross-Correlation, the problem is when recovering the text which is highly affected.

As I said before the text is binarized so when I perform an attack it is very likely that some bits representing the text are flipped. If this were a binary image and not text, the extraction would be almost imperceptible to human eye but when dealing with text everything screws up.

So my question are: do you there another way to embed text different from binarization? Do you know any other algorithm to hide text on audio that resists common audio attacks? Any other recommendation to hide same data on audio file?

Thank you


texts: Question on the update step of Bayesian filter

texts:
The update step in Bayesian filtering is given as follows, in a textbook

$p(x_k|y_{1:k}) = \frac{1}{Z_k}p(y_k | x_k)p(x_k|y_{1:k-1})$,

$Z_k$ is defined as $\int p(y_k|x_{k})p(x_k|y_{1:k-1})dx_k $

I was able to derive the second term on the right .The steps are given below
$p(x_k|y_{1:k}) = \frac{p(x_k , y_{1:k})}{p( y_{1:k})}=\frac{p(y_{1:k}|x_k)p(x_k)}{p( y_{1:k})} $

$=\frac{p(y_k,y_{1:k-1}|x_k)p(x_k)}{p( y_{1:k})}$

$=\frac{p(y_k|y_{1:k-1},x_k)p(y_{1:k-1}|x_k)p(x_k)}{p( y_{1:k})}$

By Markov property and Bayes rule

$=\frac{p(y_k|x_k)p(x_k|y_{1:k-1})p(y_{1:k-1})}{p( y_{1:k})}$
$=\frac{p(y_k|x_k)p(x_k|y_{1:k-1})p(y_{1:k-1})}{p(y_k, y_{1:k-1})}$

Applying bayes rule to the denominator

$=\frac{p(y_k|x_k)p(x_k|y_{1:k-1})p(y_{1:k-1})}{p(y_k| y_{1:k-1})p(y_{1:k-1})}$

$=\frac{p(y_k|x_k)p(x_k|y_{1:k-1})}{p(y_k| y_{1:k-1})}$

Can someone explain how the denominator in the above equation can be written as

$Z_k$ = $\int p(y_k|x_{k})p(x_k|y_{1:k-1})dx_k $


texts: Help with Kalman Filter implementation for estimating 3D position

texts:
I wrote a kalman Filter implementation using the Eigen Library in C++ and also using the implementation at this link to test my filter: My prediction step looks like this:

void KalmanFilter::Predict()
{
    // state Estimate = state transition matrix * previous state
    // No control input present.
    x = A * x;

   // State Covariance Matrix = (State Transition Matrix * Previous State
   Covariance matrix * (State Transition Matrix)^T ) + Process Noise

   P = A * P * A.transpose() + Q;
}

while my update step is:

void KalmanFilter::Update(VectorXd z)
{
    //Kalman Gain = (State Covariance Matrix * Measurement matrix.transpose) * (H*P*H^T + Measurement Noise)^-1
    K = (P * H.transpose()) * (H * P * H.transpose()+ R).inverse();

    //Estimated Stated = Estimated state + Kalman Gain (Measurement Innovation)
    x = x + K*(z - H * x);

   //State Covariance matrix = (Identity Matrix of the size of  x.size * x.size) - K* H * P;
    long x_size = x.size();
    MatrixXd I = MatrixXd::Identity(x_size, x_size);
    P = (I - K * H) * P ;

}

My initial values are:

pos_x = 0.0;
pos_y = 0.0;
pos_z = 1.0;
vel_x = 10.0;
vel_y = 0.0;
vel_z = 0.0;
acc_x = 0.0;
acc_y = 0.0;
acc_z = -9.81;

and I'm generating "fake data" by doing the following in a loop:

double c = 0.1; // Drag resistance coefficient
double damping = 0.9 ; // Damping

double sigma_position = 0.1 ; // position_noise

// Create simulated position data
for (int i = 0; i < N; i ++)
{
    acc_x = -c * pow(vel_x, 2);  // calculate acceleration ( Drag Resistance)

    vel_x += acc_x * dt;  // Integrate acceleration to give you velocity in the x axis.

    pos_x += vel_x * dt; // Integrate velocity to return the position in the x axis

    acc_z = -9.806 + c * pow(vel_z, 2); // Gravitation + Drag

    vel_z += acc_z * dt;   // z axis velocity

    pos_z += vel_z * dt;  // position in z axis

    // generate y position here later.

    if(pos_z < 0.01)
    {
        vel_z = -vel_z * damping;
        pos_z += vel_z * dt;
    }

    if (vel_x < 0.1)
    {
        acc_x = 0.0;
        acc_z = 0.0;
    }

    // add some noise
    pos_x = pos_x + sigma_position * process_noise(generator);
    pos_y = pos_y + sigma_position * process_noise(generator);
    pos_z = pos_z + sigma_position * process_noise(generator);

I then run my prediction and update step by:

// Prediction Step
kalmanFilter.Predict();

// Correction Step
kalmanFilter.Update(z);

where z is a 3 x 1 vector containing pos_x, pos_y and pos_z

My State Transition Matrix A looks like this:

and P is:

R is a 3 x 3 Matrix of:

and Q is G * G.transpose()* a * a;
where G is a 9 x 1 Matrix of

 and

a is 0.1 ( acceleration process noise)

My issue is my estimated position for y and z are off and diverge from the "real" positions. If you look at the following graphs,

This is what pos_x looks like:

This is what pos_y looks like:

And finally Z:

This is my first foray with Kalman filters and I'm not sure what I'm doing wrong here. My final goal would be to use this to estimate the position of a drone. Additionally, I have the following questions:

In a real life situation for a drone for example, how do you about choosing your Process Noise if you can't directly observe the process? Do you simply just pick arbitrary values?

My apologies for the long post. Any help is appreciated.


texts: What's the difference between Channel shortening and time-domain equalizer

texts:
I'm asking about the difference about channel shortening and time-domain equalization. I have checked on-line and tried to get something clear but I couldn't.

What I found is that channel shortening means almost  "Calculating the coefficients for DFE equalizer" But I couldn't also imagine that.

Could you please explain.

Thank you


texts: Reduce the "glare" effect of a receiver within a system

texts:
I'm working this summer on a very interesting, albeit difficult to understand at 100 %, project and it turns around reducing the glare effect than could be applied on a receiver when he gets signals coming from two emitters sending $X(t)$ (more energy for this signal compared to the other) and $s(t)$ (much weaker).

Let me put it in other words.

The receiver gets a signal mixture $Y(t)=H(t)*X(t)+g(t)*s(t)$ but we don't know $H$ nor $X$ nor $g$ nor $s$ nor the frequencies of the emitted signals. The only thing we know is $Y(t)$ and the fact that their spectra are dissociated in the best cases.

A picture speaks a thousand words, so...

Basically, the signal $X(t)$ overpowers $s(t)$. So I was thinking of a way to get the $s(t)$ and the $X(t)$ and with the few info I've gotten, I couldn't implement the generic solutions. However, I still have one card which is the FASTICA but that would suppose having many samples of "mixed " signal to extract the two important signals. What do you think about it?  I'm open to any software and/or mechanic suggestion.

I've read an interesting paper on this but they were talking about the near-far effect for GPS signals.
Thanks for the help!!


texts: Independence of Noise at Each DFT Output

texts:
My math may be a little rusty, so I would like confirmation or correction or my calculations here.

Given white noise samples, $x_i$, which are IID and zero-mean, and variance $\sigma^2_x$. I want to show that the outputs of any pair of DFT filters is also IID.

Knowing covariance of the outputs should be enough to draw a conclusion. So I examine the following:

Given the outputs $X_k$ and $X_{k'}$ defined as

$$X_k=\Sigma_{n=0}^{N-1}x_ne^{-2\pi ikn/N}$$

$$X_{k'}=\Sigma_{m=0}^{N-1}x_me^{-2\pi i{k'}m/N}$$

We can calculate the expected value of $X_k X^*_{k'}$

$$\begin{eqnarray*}
E(X_k X^*_{k'})
&=&E\left(\sum_{n=0}^{N-1}x_ne^{-2\pi ikn/N}\sum_{m=0}^{N-1}x_me^{2\pi i{k'}m/N}\right)\\
&=&E\left(\sum_{n=0}^{N-1}\sum_{m=0}^{N-1}x_n x^*_m e^{-2\pi i(kn-{k'}m)/N}\right)\\
&=&\sum_{n=0}^{N-1}\sum_{m=0}^{N-1}E\left(x_n x^*_m \right) e^{-2\pi i(kn-{k'}m)/N}\\
\end{eqnarray*}$$

We know that $E(x_k x^*_{k'})=\sigma^2_x$ when $k={k'}$ and $E(x_k x^*_{k'})=0$ when $k\neq {k'}$. So we can simplify:

$$\begin{eqnarray*}
E(X_k X^*_{k'})
&=&\sigma^2_x\sum_{n=0}^{N-1}e^{-2\pi i(k-{k'})n/N}\\
&=&N\sigma^2_x|_{k={k'}}
\end{eqnarray*}$$

But if I'm right so far, I still am fuzzy on the case where $k \neq {k'}$. I expect it to equate to zero. Can someone help me out with that part?


texts: (Newbie question) Calculation DC-offset and RMS

texts:
Being a newbie on this topic I beg you to be indulgent to me
if my questions should be kind of "too" simple ;-)
For doing automatic analysis on my audio material
I started calculating the most simple values from the int16 samples:
Clipping count, DC-offset and RMS.
For manual manipulating audio files I use an older version
of Audition (v.3), ... but it returns different values
from the ones I calculate.
I tried Audacity - but it causes my soundcard's relais to
produce strange sounds when starting - so I better leave that.
Also tried sonic visualizer (with bbc analyse plugin),
but the latter crashes the app.
So my question:

Is there a free app for analysing audio
(amplitude, rms, peaks, DC-offset, frequencies etc.)?

If you have a look at the audio piece.
http://www.mediafire.com/file/4bgf9yf6myxmfg3/test_mono.flac/file
Calculating DC-offset is the most trivial
(average of the samples), returned as % of 32768.
My result: 0,12. Audition: 0,059.

Who is right?

Clipping count
For that I use a kind of heuristics (derivation of samples),
which gives me similar results as Audition.
RMS
Audition returns different "aspects" of that.
(I don't have the English version, so I translate it myself)
Overall:

I calculate that as average of all samples' RMS, correct?
(My value: -8,6 dB. Audition: -8,97 dB)

The other RMS aspects are dependent of the window size.
I do understand it as "parts" of the whole audio,
in which the RMS is calculated indepently.
Taking all parts' values you can take
the min, max and average values of them.

Correct?

But strange enough, Audition's help description of that
dialog advices to take a big window size for high dynamic range
audio (and vice versa).

I would have thought the other way round (???)

OK (using 50ms for window size).
Min (my app: -23,1 dB. Audition: -inf).
Max (my app: -4,8 dB. Audition: -4,5 dB).
Average (my app: -9.4 dB. Audition: -9,9 dB).

Who's wrong?

Last but not least I "invented" a criterion for dynamic range
(in sense of "loudness war", dynamic compressed audio),
and calculated the standard deviation of all the windows' RMS.
As percentage of the average RMS.
(Gives 43,4% for my sample)

Does this make sense, for a simple approach?

I hope my questions don't hurt too much ;-)
And maybe answered by few short hints.
Thanks in advance.


texts:
I have written a small simulator where I am trying to simulate the MRI signal generation process. I have a SINC RF Pulse to excite the spins and a readout gradient and k-space sampling during the readout. The k-space is sampled at equidistant points in time and the signal is inverse fourier transformed to give the image.

The synthetic image phantom for which I am trying to generate the signal is completely homogeneous. So, I have spin objects which have identical T1 and T2 relaxation times and density.

So, when I do the reconstruction, I was expecting the same intensity at each of these sample points. Instead what I get is something as follows:

Now, I am not sure if this is because of an error I made or simply because of numerical issues or the fact that the SINC RF pulse is terminated. I also changed the RF pulse to simulate a perfect flip i.e. all spins end up along the Y-axes perfectly, but this effect is still seen.

I wonder if this is due to Gibbs artefacts? I have a feeling this has something to do with how the fourier transform works but I am not sure.

I increased the number of k-space samples and also increased the FOV to 0.5 meters and then I get the following plot:

Finally, I did a reconstruction with 128 sampling points and a FOV of 0.128 meters and I got the following image:

P.S: Do not worry about the y-values as I was plotting these with different spin density values.

[EDIT]

Here is some more detail on the simulation. The underlying imaging object constitutes a single k-space line for simplicity. All the spins have the same T1, T2 and spin density with values 0.833, 0.083 and 1 respectively.

The FOV of the underlying object os 0.128 m. The constitutes a k-space matrix of 128 elements, 1 mm apart.

The RF excitation was done for 2 ms and the centre frequency was set to 0 Hz and a bandwidth of 2000 Hz. The flip angle was 90 degrees.

Immediately after the RF pulse a rewinder gradient was applied to dephase the spins. This was done to ensure the echo happens in the centre of the k-space line. The rewinder gradient had the peak amplitude of -0.022 T/m. The duration of this was 0.000816 sec and the area under the gradient was about -1.3934e-05.

Immediately after the rewinder, the readout gradient was applied with a duration of 0.0014 seconds and the area under the gradient was 2.7868e-05. The area of this gradient was twice that of the rewinder gradient to allow for the echo to form in the middle.

During the readout gradient, the k-space was sampled at uniform intervals. The transverse magnetisation over all spins was summed together at each sampling time point to yield one complex k-space value.

This was inverse fourier transformed using the scipy python package as:

np.abs(fftpack.ifftshift(fftpack.ifftn(signal)))

So, it shows the absolute value of the signal.

Also, more specifically, there is no spoiling performed as only one k-space line is read and only one RF-excitation is performed.

The simulation is performed by solving the Bloch equations using a fourth order Runge-Kutta method with a step size of 10e-6 sec. I see the same behaviour when I decreased the step size.

So, once I perform the simulation, I detect an echo as follows:

and the reconstructed image is the plot here:

Note that because the object is only simulated in one-dimensional, the image is just a 1-D plot.

My issue is I was expecting a k-space line which was close to a spike in the middle since the object is completely uniform.

[EDIT 2]

I figured out that I had made a small mistake in the readout part (I was not sampling at the correct time). Fixing this makes the echo less zagged but it is still not perfect. The echo and the reconstruction is below:


texts: Recursive systematic convolutional code (RSC) realization in MATLAB

texts:
How can I implement $\rm RSC(2,1,2)$ in MATLAB? Can I use convenc function to build RSC?

The convenc function has 2 arguments: code = convenc(msg,trellis) where trellis=poly2trellis(ConstraintLength,CodeGenerator), and the generator polynomial for this RSC is:
$$
G=\left(1,\frac{1+D^{2}}{1+D+D^{3}}\right)\quad\text{or}\quad G=\left(1,\frac{5}{7}\right)
$$
but for poly2trellis argument CodeGenerator need to be in octal.


texts: disadvantages of FFT, it can not extract enough frequencies without enough samples

texts:
Let's say sampling rate is $Fs = 44\mathtt{kHz}$, now I have $N = 2048$ samples, then I can get $N/2 + 1 = 1025$ frequencies.

I'm confused by Matlab's FFT documentation that says the frequencies are present at:

(Fs/2) * linspace(0, 1, N/2+1)

It is very probable that the real frequencies are not among  the frequencies generated. What should we do to get the real frequencies?


texts: Feature Detection on a Circular Metal Surface for Character Recognition

texts:
Hello I'm trying to recognize characters on the circular metal surface.

Without any preprocessing I could match cropped features with  original image.
I've actually cropped 2 characters together to increase interested points on a template.

however on different locations of the conveyor belt and with different illumination angle and perspective SURF matching fails.

I will definately crop the ROI from conveyor belt but still I am not sure whether I have to create a database of images to increase accuracy or some kind of preprocessing (I've read that for feature detection there is no require preprocessing on images)

I am open to new suggestions for this problem. Maybe a different algorithm or method that can help me out.

thanks in advance.


texts: How to install FCLAB for EEG connectivity calculation

texts:
I am doing about another person question here, so the question is:

How I could do EEG functional connectivity calculation by Matlab?

So I have googled the EEG functionality and find two plugin by the name of SIFT and FCLAB.
Also, I checked FCLAB installing by google and find out one instruction in the following link:
https://github.com/ramsys28/FCLAB/issues/1
Which said copy FCLAB folder to the EEGLAB plugin folder, FCLAB files can also be found at the following link:
https://github.com/ramsys28/FCLAB
BUT my OS is Linux base which I not have Matlab on it, and I working more with python software codes, so I asked here about that the above instruction is working or not?
Thanks for your attention.


texts: How to change the angle unit labels on polar plot in MATLAB?

texts: I hope this is not off-topic. I'm wondering is there anyway to change the default angle unit label in MATLAB polar plot. The default angle label is from 0 to 360, at 30 increment. It goes like 0,30,60,... But I want to change it, say, from 0 to -360, at -30 increment. So the labels will be like 0,-30,-60,... Is there anyway to do it?


texts: Calculating the energy of trapezoid wave

texts:
Here is the trapezoid wave,although i know we can use $\int^{4}_0 |x(t)|^2dt$ to calculate its power,but i want to divide this trapezoid to right triangle$(t=0$~$2)$,and a square ($(t=2$~$4)$).

I know the the triangle effective voltage is $V_{eff}=\frac{V}{\sqrt{3}}$,and i want to use $P=\frac{V_{eff}^2}{R}$ to calculate the power ,  but i am not sure about which  triangle is correct,i mean,for which triangle(a or b or c),i can use $\frac{V}{\sqrt{3}}$ to calculate its power,or are these three wave's power the same,but different in energy?


texts:
I want to study stress based on sound recordings. However, there are many concerns about keeping the identity of the recorded persons anonymous.
I want to keep intensity and spectral contents intact but I am not interested in the message or the identity of speakers.
I thought about time reversing the recordings, but this could be easily undone by someone.
Do you have any suggestions for this purpose?


texts: What does mean and variance do in Gaussian noise?

texts:
I = imread('us-currency.png');
mean = 0.5;
variance = 0.5;
J = imnoise(I,'gaussian', mean, variance);
imshow(J);

What are the effects of mean and variance on that image?

Does mean increases the brightness?
Does variance increases the amount of noise?

Input

Output


texts: Why are the Continuous Wavelet Transforms of the same signal drastically different?

texts:
I'm currently studying wavelets and am running into confusion with regards to CWT coefficients. Ideally, I want a CWT algorithm that produces outputs similar to that of a STFT - i.e. produces coefficients that correspond to frequency magnitudes of the input signal. However, I'm finding conflicting results.
From the CWT Wikipedia page, I'm given the illustration below. This illustration shows that the wavelet coefficients are essentially just sampling different subbands of the signal. This is not what I'm looking for.

However, when I look at a MATLAB web page discussing the CWT and DWT, I find a different illustration of the CWT when operating on the same input signal. This diagram looks more similar to what I was hoping for where each frequency has a constant coefficient.

What is the difference between CWTs?


texts: Audio wave amplitude decreased after FIR filter

texts:
I am doing sampling and filtering for a 17KHz sine audio wave signal, the Nyquist frequency is about 20.8889KHz. I found the fact that this audio signal's amplitude decreased dramatically after it went through an FIR bandpass filter( pass band 16.8khz-17.2khz), the amplitude also decreased after an FIR highpass filter(cutoff frequency 17khz). I think that's because the audio signal frequency is kinda close to Nyquist frequency, am I right?
Is there a filter can preserve the original amplitude of the audio signal within the frequency I want?
Shown below is a pic consists of 2 signals before and after filtering(highpass filter)
Thank you in advance


texts: Is step detection the correct approach to this problem? what if not?

texts:
I'm looking for some advice on where and what to start reading for learning to solve this.

I've the time series of the position coordinates (x,y) of an animal in an open field (just a cage). I want to detect the time instants when the animal stops or starts "walking" (i.e.: moving from one place to another, not in the same place).
I  thought it couldnt be so difficult, but I'm finding some trouble solving it.

So, I'm looking for the times from which the signal will be flat, or the instant where the signal starts changing from flatness:

I think this thing may be related with step detection but I'm not sure.
I would start with just one coordinate for simplicity. Is seems that step detection could be what im looking for, but:
a. my signal "baseline" will be different every time. the animal will move and stop somewhere else, and so.
b. the signal can be VERY "noisy" since the animal could and will move a lot while staying in the same place (ie: grooming)
c. this changes can be whether slow or fast, and I need both.

Firstly, i ll be glad if i can solve this problem for just one coordinate, although I will have to look for this changes in the two coordinate system.

so, my questions are:
1. is step detection a good aproach to this problem? what else if not?
2. any suggestion in doing this for both (x,y) coordinates?

thanks in advance

EDIT: I get (x,y) coordinates by acquiring an overhead image with a camera and tracking a led attached to the animal's head. Tracking is done by color filtering.

EDIT2: copy of the data:

https://www.dropbox.com/s/oph33szu891rgrl/data_sample.txt

https://www.dropbox.com/s/idulfdr965eeh7i/data_sample.mat

data format is (x,y,t)

EDIT3:
I've been trying smoothing the data but it is not really what I expected. I need to clean not high frequencies but low amplitudes. The movements I look for can be either fast or slow, but with big amplitude.
here, an example of signal with matlab function smooth()
note that I'm looking for the green moments, which I lose when smoothing


texts:
How can I prove the equation $(3)$? I can't understand why there is a $2/N$ in $(3)$. Why he just get the second term in DFT?

  Consider a sinusoidal input signal of frequency $\omega$ given by

  $$
x(t) = \sqrt 2 X \sin\left(\omega t + \phi\right)\tag{1}
$$

  This signal is conventionally represented by a phasor (a complex number) $\bar X$

  $$
\bar X = Xe^{j\phi}=X\cos\phi+jX\sin\phi\tag{2}
$$

  Assuming that $x(t)$ is sampled $N$ times per cycle of the $60\textrm{ Hz}$ waveform to produce the sampled set $\left\{x_k\right\}$

  $$
x_k = \sqrt 2 X\sin\left(\frac{2\pi}{N}k + \phi\right)\tag{3}
$$

  The Discrete Fourier Transform of $\left\{x_k\right\}$ contains a fundamental frequency component given by

  $$
\color{red}{\boxed{\color{black}{\bar X_1= \frac 2N \sum_{k=0}^N x_k e^{-j\frac{2\pi}{N}k}}}}\tag{4}\\
$$


texts: Mismatched offset when computing cross-correlation using FFT

texts:
In MATLAB if I compute the cross-correlation of 2 arbitrary signals using the xcorr function I get exactly the answer I expect. If I then calculate the cross-correlation using element-wise fft multiplication, the resulting signal is incorrect (offset) by 1 element. Can anyone tell me what I'm missing? Here is my MATLAB script:

a = [0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0];
b = [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0];
elems = length(a);

% Calculate the cross-correlation coefficients using built-in function
[ac_coeffs, lags] = xcorr(a, b);
subplot(2, 1, 1);
plot(lags, ac_coeffs);

% Now do it manually by fft multiplication
a = [a, zeros(1, elems - 1)];
b = [b, zeros(1, elems - 1)];
af = fft(a);
bf = fft(b);
mult = af.* conj(bf);
result = ifft(mult);
subplot(2,1,2);
plot(result);


texts:
I wanted to do formant analysis using MATLAB and did everything according to the following tutorial:
https://de.mathworks.com/help/signal/ug/formant-estimation-with-lpc-coefficients.html
I have speech signals with a sampling frequency of 44.1kHz, therefore I changed the order of the LPC estimation to 46.
But I have the following problem: the algorithm sometimes works well, but most of the time the formants are completely off. Especially for female voice. I tried to play around with the required bandwidth and lower limit for the formant frequency, but nothing yielded sufficiently good results for various voices.
I hope, that someone might have some tips for me to improve the code, such that it works more generally.


texts: Is Smith's book still usable and up-to-date?

texts: I'm learning DSP with Smith's book, but I discovered the second edition is from 1999, is this book still usable today?


texts: Calculation of gain at certain frequency in case of nth order IIR filter


texts: Discrete-time Convolution Convergence Issue

texts:
I'm trying to solve a problem on convolution from Alan V.Oppenheim:
Find the convolution output $y[n]$ for the following signals:
$$x[n]= u[n]\quad\text{and}\quad h[n]=a^{-n}u[-n-1], \ a>1 $$
I started the evaluation:
$$y[n]=\sum_{k=-\infty}^{+\infty} u[k]a^{-n+k}u[-n+k-1]$$
considering that $u[k]=1$ for $k\ge0$ and $u[-n+k-1]=1$ for $k\ge n+1$
which I evaluated to
$$y[n]=a^{-n}\sum_{k=n+1}^{+\infty} a^k$$ where $n$ could be $<0$ or $>0$
and I tried to evaluate for $n>0$:
where I faced an issue:
I subsittuted $k-1$ to $m$ and reframed the equation as:
$$y[n]=a^{-n+1}\sum_{m=n}^{+\infty} a^m=a^{-n+1}\left[\left(\sum_{m=0}^{+\infty} a^m\right)-\left(\sum_{m=0}^{n-1} a^m\right)\right]$$
now for $a>1, \ \displaystyle \sum_{m=0}^{+\infty} a^m$ will not converge. How will I evaluate this?


texts:
I have a .wav file which came off of a cd with a scratch, so there is a periodic "clicking" of about 6 Hz.  I was playing around with trying to squash the clicking by taking FFT and filtering the range around this frequency, but with no success.  I tried smoothing the low end of the FFT spectrum, but it didn't change anything.  Then I tried zeroing this part of the spectrum up to higher and higher indices, but it still wouldn't kill the click.  Eventually it started to attenuate the entire signal, but the click was still there, and perhaps even more prominent than before (relative to the rest of the signal).  Why is this the case?

My first speculation is that the click is somehow not quite steady enough to be localized in a single region of the spectrum.  Can anyone confirm this hypothesis or else explain the true source of the issue?  Also, is there a better means of eliminating the click?


texts: Deconvolution of an Image Acquired by a Square Uniform Detector

texts:
So, I acquired some images by scanning a radiation source with a square detector like in the following gif.

Where the dashed grid represents reality, the 3x3 square my detector, and the 4x4 my acquired data.
Obviously, because I'm moving the detector in steps smaller than its size, my image is convoluted. Assuming my detector has uniform sensitivity, how do I go about deconvolving it?

The point spread function should be simple, but I can't seem to find examples anywhere.

(I can work with Python or Matlab, so feel free to give any examples you want)

The data I'm working with consists of gaussian-ish peaks surrounded by noise like so:


texts: Fourier transform as the integral of a parameter multiplied by an homogeneous wave

texts:
Can a Fourier transform in space be interpreted as the integral of a parameter multiplied by an homogeneous wave $\sigma$?
where $\sigma$ is:
$\sigma$=$e^{-ikx}$
Are there papers or book that illustrates this interpretation?


texts:
I would like to down sample a signal from it's original sampling rate of 2,000,000 samples per second to 48,000 (or 44100) starting from the output of an FFT.  This FFT delivers 16384 data points (8192 real and 8192 imaginary).  The data comes from a software defined radio as I and Q data.
I cannot do this in the time domain as I do not have the tools in my app to do so.  So I need to start from the FFT, drop the necessary data, then perform an iFFT.  These tools I have.
What I need help with, is determining how many samples to drop and how to arrange the data afterwards to be delivered to an inverse FFT which needs real and imaginary data.  I understand that you just drop the upper frequencies.  But then what?
Can someone help out please?
Thanks