Why does Hilbert filter distorts the shape of the signal? If all the harmonics composing a the signal are shifted by the same amount, this would be the same as sampling later or earlier in time. I think.

Take a simple pulse train as an example.
If Fourier is performed when the pulse is high for time t after 0, and flat for time t before the 0 it will have odd symmetry and only sin harmonics. If the same signal was sampled differently with half of the pulse left of y-axis and other half right of it, it would have even symmetry and contain only cos terms.

Each time actual signal is the same, just its Fourier representation is different.

Difference between these is only shifting all harmonics by 90 degrees. Why is Hilbert filter different in those regards?

I know square wave passed through Hilbert filter produces sideways curly braces like shape, so my fundamental understanding of something is wrong. Thanks for help.

A Hilbert transform is dispersive, as in any constant phase shift over frequency, which will distort a signal that has bandwidth. This is because a fixed DELAY, which would have no distortion, has linear phase versus frequency. If you fix the phase versus frequency to be constant (as we do in a Hilbert transform), the delay at each frequency will be different and thus will cause significant distortion to the signal.

This is demonstrated with my slide below that is used for explaining why linear phase filters are important specifically, but also helps demonstrate what a constant phase shift versus frequency would do. Note that delay specifically is the derivative of phase versus frequency:

$$delay = \frac{d\phi}{d\omega}$$

In the slide below, the first three Fourier components of a square wave signal are shown, along with their addition when they all add up in time (which would happen with a linear phase filter as the delay in time given the formula above would be the same for each frequency). The "Non-linear Phase" case shows what happens if the Fourier components are summed after going through different time delays, which occurs in non-linear phase filters, AND occurs if the composite signal goes through a system where the phase is constant for all frequencies (as the time delay for each frequency must be different in order for the phase to be the same).


Create a third octave spectrum from a time signal I have to create a third-octave spectrum from a time signal on Octave GNU. I found some code on the net to help me but I don't have all the parts of the algorithm.

I have a .csv file which contains a simple sinus temporal signal. The specifications are:

Duration: 0 to 180s
Frequency: 32768

I calculate the octave bands from 10Hz to 10kHz :

    fMin=10;
    fMax=10000;
    octs=log2(fMax/fMin);
    bmax=ceil(octs/bw);

    %Octave bands.
    fc=fMin*2.^((0:bmax)*bw);    %Center frequencies.
    fl=fc*2^(-bw/2);             %Lower frequencies.
    fu=fc*2^(+bw/2);             %Upper frequencies

After I can display all the third octave filters using Butterworth filters :

    bw = 1/3;
    numBands = length(fc);

    b = cell(numBands,1);
    a = cell(numBands,1);

    figure
    for nn=1:length(fc)
      [b{nn}, a{nn}] = butter(2,[fl(nn) fu(nn)]/(fs/2));
      [h,f] = freqz(b{nn}, a{nn}, 1024, fs);

      hold on;
      plot(f, 20*log10(abs(h)));
    end
    set(gca, 'XScale', 'log')
    ylim([-50 0])

Which gives to me :

The last filters are great but the first ones have missing points.
I don't know how to fix that and, mainly, how can I use this to make my third octave spectrum.
Does someone knows how to do this ?

Hate to answer this question, but it is a plotting issue, not the actual design problem. Basically, you don't have enough points at low frequencies cause they are spaced linearly and you are plotting on logarithmic scale. You could do:

[h,f] = freqz(b{nn}, a{nn}, logspace(log10(fMin), log10(fMax), 1000), fs);

Explicitly specifying the logarithmically spaced bins. Or you can use a very large N (something like $2^{16}$).


How to solve a two variable equation of images? I have two equations of images

$$g_1=f_1+h_2*f_2$$

and

$$g_2=f_2+h_1*f_1$$

where $*$ is convolution and $g_1,g_2,h_1,h_2$ are known.

How do I find $f_1$ and $f_2$?

If we add hats to denote the Fourier transforms of those functions, convolution turns into multiplication, and addition stays as addition:

$$\hat g_1=\hat f_1+\hat h_2\times\hat f_2$$
$$\hat g_2=\hat f_2+\hat h_1\times\hat f_1$$

If we feed that to Wolfram alpha, sans hats:

solve({g_1=f_1+h_2*f_2, g_2=f_2+h_1*f_1}, {f_1, f_2})

the output is, with hats added back:

$$\hat f_1 = \frac{\hat g_2 \hat h_2-\hat g_1}{\hat h_1 \hat h_2-1} \text{ and } \hat f_2 = \frac{\hat g_1 \hat h_1-\hat g_2}{\hat h_1 \hat h_2-1} \text{ and } \hat h_1 \hat h_2-1\ne0$$

Frequency domain division equals time domain deconvolution. The frequency domain $1$ represents a time domain identity convolution kernel of a single 1 at $(0, 0)$ and all other values zero. Reading from the above, the solution is valid only if $h_2$ convolved by $h_2$ minus the identity kernel gives an invertible filter, in other words one without frequency domain zeros.


What happens to the phase spectrum when I resample? [An update is added at the end of the post after receiving first response]
I have an algorithm which is very sensitive to phase shifts. It works with signal sampled at 40MHz (it's a neural network so actually the sampling rate is not very flexible)
I have a data sampled at 21.818e6 MHz and I resample it to ~40 MHz by:
sample_cnt = len(test_data)
Fs_org = 21.818e6
Fr = 40e6
resam_cnt = np.floor(Fr * sample_cnt / Fs_org)
test_data_re = scipy.signal.resample(test_data, resam_cnt.astype(int), axis=0)

test_data is sensor data with 1704 sample numbers so (1704,1).
When I look at frequency spectrum and phase response of the resampled data the frequency spectrum matches original signal (it's a bandpass signal with center frequency of 5MHz) but phase response for higher frequencies is flat (see image).

pyolt.figure()
pyplt.phase_spectrum(test_data_re, Fr, fc=5e6)
pyplt.phase_spectrum(test_data, Fs_org, fc=5e6)

Blue line is corresponded to the phase response of resampled data and orange line is corresponded to original data which is overlapped on the resampled data.
1- Is this resampling correct?
2- Is there a solution to create any kind of corresponding phase information for higher frequencies?
[Update]
I resampled same data with Matlab:
resample(test_data,2,1)

Then saved the data, loaded in python and compared two resampled signal (see image)

The blue one is the python resampled data and the orange one is the matlab resampled data
I do agree with @Hilmar that "When you re-sample you extend the frequency range, but the energy at these new frequencies (between 11 MHz and 20 MHz) is close to or exactly zero. That's the whole point of a "good" interpolation: make sure you don't add stuff that wasn't there before.
That means the phase at the new frequencies is either undefined or doesn't really matter since there is no relevant energy. If the magnitude is zero, it doesn't matter what the phase is"
Now it rises this question that how come Matlab is creating phase information in those frequencies?

1- Is this resampling correct?

No. Resampling is almost never 100% correct but it usually requires a complicated trade off between different artifacts and properties: passband extension and flatness, stop-band attenuation and size, aliasing rejection, phase distortion, onset & transient preservation, causality, CPU and memory requirements etc.
Which trade off is the best, depends a lot on the requirements of your specific application.

2- Is there a solution to create any kind of corresponding phase information for higher frequencies?

Why would you need that? Your original signal cannot contain any energy at frequencies equal or above the original Nyquist frequency. When you re-sample you extend the frequency range, but the energy at these new frequencies (between 11 MHz and 20 MHz) is close to or exactly zero. That's the whole point of a "good" interpolation: make sure you don't add stuff that wasn't there before.
That means the phase at the new frequencies is either undefined or doesn't really matter since there is no relevant energy. If the magnitude is zero, it doesn't matter what the phase is.


Finding impulse response of an RC circuit from its step response I have to find the impulse response of an RC circuit (c up ). I have to find it from the step response $g(t)$.

I found that

$V_i - V_c - V_r = 0$ , with $V_r = V_u$

So

$V_u = V_r = R I_r$ , with $I_r = I_c$

$V_u= RC  \frac{dV_c}{dt}$

$V_u = y(t) = RC \frac{d(x(t)-y(t))}{ dt}$

because $V_c = x(t)-y(t)$

From this i found $y(t) = RC ( x'(t) - '(t) )$ , i simply rewrote d/dt with '.

Now if i want to found the step response, i considered $x(t) = u(t)$ and $y(t) = g(t)$ that is the step response

So i wrote
$$g(t) = RC u'(t) - RC g'(t)$$

but $u'(t) = 0$ so is $RC g’(t) + g(t) = 0$

Solving this differential equation i obtained that

$$g(t) = k_0 +  k_1 \cdot e^{-t/RC }$$

Until now i obtained the same results of my professor but now that i have the step response i have to find the $h(t)$, impulsive response , but I don’t know how. Can someone please help me ? This is the same problem of my previous question but in this case i have to find the impulsive response using the step response. Thank you so much

The answer (How to calculate the impulse response of an RC circuit using time-domain method)  provides a direct time-domain solution of an RC circuit for the impulse reponse $h(t)$. Now this new answer modifies it to solve for the step-response $s(t)$ instead and then computes the impulse response according to :

$$h(t) = s(t)' $$

The differential equation of the first order circuit was derived as
$$ \boxed{ y'(t) +\frac{1}{RC} y = x'(t) } \tag{1} $$

The step reponse $s(t)$ is defined as the output $y(t)$ of Eq.(1) when the input $x(t)$ is a unit-step function $$ x(t) = u(t) \implies y(t) = s(t)  $$

Let's apply a one stage direct solution to obtain $s(t)$.

The homogeneous solution is found from
$$ y'(t) +\frac{1}{RC} y = 0  \tag{2} $$

The characteristic equation : $s + \frac{1}{RC} = 0 \implies s = - \frac{1}{RC}$.

The (causal) homogeneous solution is :

$$y_h(t) = K e^{-t/RC} u(t) \tag{3}$$

Then, the particular solution $y_p(t)$ will be from the method of undetermined coefficients as follows:

For the particular input $x(t) = u(t)$ we may assume a particular solution as $y_p(t) = A u(t) + B \delta(t)$, then plug this assumed solution into Eq.(1) to find out the coefficients $A$ and $B$. This yields $A=0$ and $B=0$, hence the particular solution is found to be identically zero.

Then since the total solution is:
$$s(t) = y_h(t) + y_p(t) = y_h(t)$$

to find the value of unknown $K$ in Eq.(4) we need an initial condition on the output $y(0^+) =s(0^+) $ ; found from the circuit physics: apply KVL around the loop of voltage source, capacitor and resistor, with the fact that the capacitor voltage is fixed at zero at time $t=0$. Then we find $y(0^+) = s(0^+) = 1$ which yields $K=1$

Hence the step-response is:
$$ \boxed{ s(t) = e^{-t/RC} u(t) } \tag{4}$$

Then the impulse response of Eq.(1) is found to be
$$ h(t) =  s(t)' = \left( e^{-t/RC} u(t) \right)' $$

$$ \boxed{ h(t) = \delta(t) - \frac{1}{RC} e^{-t/RC} u(t) } \tag{5}$$


Testing if a random signal was modified or if it contains noise Given a random signal, is there a method of finding if it was modified or if it contains noise if you don't know the original signal?
E.g.: given a random image, can I find out if it was modified if I don't know the original?

No.  You have to start with some prior knowledge of the signal and of the noise.  Your signal may be my noise, and visa-versa.


zero-phasing frequency components while keeping the same magnitude, in Matlab How is it possible?

I was thinking of taking just the real part of the DFT of my signal (isn't that zeroing out phases?) with

real(fft(X))

but the magnitudes of that result don't match :

abs(real(fft(X))) != abs(fft(X))

Thanks in advance.

it is obvious that abs(real(fft(x)))!=abs(fft(x)) since you are removing the imaginary part while taking the real function.One way if you want a signal with all dft components being real.Is to normal to just take abs(fft(x)) as the spectrum itself rather than taking the real components.In this you maintain the magnitude and all the components have zero phase.I hope this answers your question


Purpose of using polyphase filter bank (PFB) What is the advantage of using a polyphase filter bank (PFB) for spectral analysis over just using the FFT? In the standard "critically sampled" uniform DFT filterbank, the polyphase decimation/filtering is followed by an $M$-point DFT block, which implements the last step of the PFB in a computationally efficient manner using the FFT.

If you have to do an FFT anyways, why bother with the PFB? Is the reason that I can choose a custom prototype low-pass filter on the front-end? Are there some computational savings I'm missing out on?

EDIT: If comparing this to a bank of quadrature downconverters, what is the point of using a PFB if the FFT is the same mathematically? It can't be the delay of $N$ samples requried to fill up an FFT block because the decimated branches have $1/N$ the rate, which means you will be waiting the same amount of time on average with either approach. What am I missing?

One advantage of a polyphase filterbank approach is, as you guessed, that you can control the frequency response of each channel. When using a DFT alone, you have limited control over the frequency band covered by each bin (characterized by a Dirichlet kernel in the unwindowed case, or by the frequency response of whichever window function you select). This is sufficient for a lot of applications.

In some other applications, however, you might want tighter control over the per-channel frequency response. Say you wanted to construct a DFT-based spectrum analyzer with very particular specifications (e.g. -3 dB response at the midpoint between output bins, -80 dB response at a one-bin spacing from center). You can utilize the polyphase filterbank structure to implement whatever filter is needed to achieve the specified level of performance.

Another application is in cases where reconstruction is needed after the channelizer: if you're careful in the design of the polyphase filter, you can actually perform straightforward spectral modification (similar to the "ideal filtering" that many signal processing beginners attempt using the DFT) in the frequency domain, then use a synthesis filterbank to take the composite signal back into the time domain. This structure, with cascaded analysis and synthesis stages, is known as a transmultiplexer.


Realtime sample rate conversion Is it possible to have a realtime sample rate conversion in a way that a peer A with an audio stream at 44.1 KHz sends it signal over the network to another peer that is on an audio stream at a 48 KHz sample rate.

If possible, how should this be accomplished? How can we workaround the fact that one peer is consuming data on a different rate than one is producing?

Thanks!

The consumption time and transmission time is identical: One second of data is still one second of data regardless of sampling time.

The greatest common divisor between the two rates is 300, thus to resample this exactly from 44.1KHz to 48KHz you would need to use the ratio $160/147$ (and the inverse for the other direction):

$147$ is factored into $3, 7^2$

$160$ is factored into $2^5, 5$

The following demonstrates one approach to resample from 44.1KHz to 48KHz, where care has been taken to not reduce the sampling rate below 44.1KHz (if that matters for fidelity concerns) and the multiple stages simplifies the filtering needed:

Interp by 4, decimate by 3, interp by 8, decimate by 7, interp by 5, decimate by 7.

This would be implemented with the following structure where the interpolator blocks signify insert of $I$ samples between each sample and the decimation blocks signify selecting every $D$th sample and throwing away the rest. The intermediate blocks can run at any arbitrary higher sampling rate to keep up with the throughput and the input/output blocks are rate matched (consuming samples at 44.1KSps and providing output samples at the 48 KSps rate).  To do this as shown with 20 MHz audio bandwidth and 80 dB resampling image rejection, I estimate that 171 taps would be needed for FIR1, 95 taps for FIR2 and 25 taps for FIR3 (as linear phase filters so one multiplier for every 2 taps).  For real-time application, the expected delay through the resampler would be 7.9 ms.

The filters could certainly be designed with windowed Sinc functions (this is known as the windowing approach to FIR filter design which is sub-optimal - see our further discussion here FIR Filter Design: Window vs Parks McClellan and Least Squares). The least squares algorithm (firls in MATLAB/Octave and Python) provides an optimal solution for resampling applications, resulting in higher image rejection for a given number of taps. Further, in many resamplers (not this one due to the close ratio), the images to be rejected can be isolated to distinct frequency bands; resulting in the use of multiband filters which the least squares algorithms support and further maximize rejection where it is needed most.

Interpolation and decimation resampling can also be accomplished by mapping the same filter coefficients as designed for the resampler shown above into polyphase structures as depicted in the diagram below. This would be identical in performance to the resampler above but can be done with an internal sampling rate as low as 67.2 KSps and significantly fewer overall computations. The decimation is done by selecting the appropriate filter output associated with the computation cycle for each decimator output rate (for the first two stages this just ends up being that the commutators move back one sample after every output update, and the last stage moves forward one sample after every other update). This can be a very efficient approach since only one of the filters in each stage actually needs to be computed for each output (note that each filter within a group would contain the exact same data but the multiply and sum only needs to be done on one of them each time). Since only one filter in each stage actually needs to be computed on any given cycle, the implementation can be done with just three FIR filters (44-tap, 12-tap and 5 tap) with the coefficients updated from a ROM table for each decimator output computation cycle (this high efficiency approach would require a tightly synchronized state machine, while instead computing all internal filters would allow for a lot of slop with the internal timing above the minimum limits). If 80 dB of image rejection or 20 MHz of audio bandwidth is not necessary, then all filter lengths can be reduced accordingly.

Resampling with polyphase filters are detailed further in this post:

How to implement Polyphase filter?

As Robert points out in the comments, the above would work continuously if the input and output were at the exact frequencies given; or if the input frequency was slightly less as the processing can ensure the next sample is ready prior to being clocked by the output. The problem will be if the input frequency error is slightly higher (or output frequency slightly lower) that lost samples will occur. Some buffering can be provided to sustain modest frequency variation but will invariably overflow. The only robust solution for a real time application without operating time limitation is to ensure the input and output clocks are synchronized through some mechanism.  Given the OP is mentioning sending the data over the network, I believe buffering would then be required based on the predicted worst case clock inaccuracies between the two locations and longest time duration of an audio transmission.  If the input and output clocks were co-located then they could be PLL locked to each other.


Filters used in FBMC On page 137, chapter 5 in the book titled 5G Physical Layer: Principles, Models and Technology, the author says:

An important attribute of the FBMC waveforms is that they typically do
not employ any CP or GI, unlike the OFDM-based waveforms. The
prototype filters in the FBMC waveforms have long decaying tails which
can make them robust to the ISI caused by frequency-selective channels

Question: Why would filter long decaying tails make a signal robust to ISI caused by frequency selectivity?

long decaying tails which can make them robust

Why would filter long decaying tails make a signal robust

can make them robust, not makes them robust!
You need a long filter length for narrow bandwidth and steep transition width.
So, with a long filter, you can make very narrow subcarrier channels.
Those being narrow, they don't "see" multipath, they only see a flat channel.


What is first order difference? I don't understand what is a first order sequence. For more context: I am reading a paper where the author extracts different sequences (in the STFT domain) and then "The resulting sequences are differentiated (via first-order difference)".

What is the meaning of first order difference here?

The first-order difference operation is a technique for numerical differentiation. It is the simplest method that I know of, and consists of just treating a discrete-time signal as piecewise linear between each pair of points and calculating the slope between them. The idea is that this time series of slopes should approximate the time derivative of the underlying signal.

$$
x_d[n] = \frac{x[n] - x[n-1]}{T}
$$

where $T$ is the sampling period of the signal $x[n]$.


Block diagram for a complex impulse response I have this question regarding digital systems which might not handle sinusoidal signals in general. So let us say that I have a system with impulse response $h(n) = \{1+j, -1, 2j \}$ with first element at index zero. Clearly this is a complex quantity. When I take the z transform,I get $$H(z) = (1+j) - z^{-1} +2jz^{-2}$$
That is $$\frac{Y(z)}{X(z)} = (1+j) - z^{-1} +2jz^{-2} \implies Y(z) = (1-z^{-1})X(z) +jX(z)[1+2z^{-2}]$$ In time domain we get $$y(n) = [x(n) - x(n-1)]+j[x(n) + 2x(n-2)]$$ So now I had to implement this as a block diagram. I will not be showing that because that is not required here. My doubt, however, is concerned with the imaginary term. What should we do with the imaginary term, this is my doubt. What I did is that I first created the signal $x(n)+2x(n-2)$ which is the imaginary part and then placed a $\frac{\pi}{2}$ phase shifter block in front of the total sum (just like a black box with a label on it). This is what we usually do in analytical signals, right. And a similar thing is done in placing in-phase and quadrature-phase components of a message signal on a carrier wave in tele-communications. I applied the same idea here. Now my professor says that this is wrong because in tele-communications we are usually dealing with sinusoidal signals, but here that might not be the case in general. What I think however is that I am still right. Even if we are in digital domain, why would it matter. A $j$ always shifts any signal by $\frac{\pi}{2}$ in phase. Am I really wrong? I do not want to prove my professor wrong but I didn't asked the reason and I am still not convinced. That is the only reason of posting this question online.

A digital system with complex input and complex output has a block diagram in which the various arrows carry complex signals, each register holds a complex number, each multiplier computes the product of two complex numbers, each adder computes the sum of two complex signals, etc. Consider, for illustration, a simple FIR filter with complex input $x$ and complex output $y$ where
$$y[n] = x[n] + a\cdot x[n-1]\tag{1}$$
with $a$ being a complex constant.  The corresponding block diagram is as shown below:

But, in real life, complex-valued signals actually consist of two separate real-valued signals traveling on two separate wires and it is only in our minds that we think these two real-valued signals as constituting a complex-valued signal (In software or MATLAB, each complex sample is stored as a pair of real numbers (two words)).  So, if we express complex input signal $x$ as $x=x_r+jx_i$ and complex output signal $y$ as $y= y_r + jy_i$ where $x_r, x_i, y_r, y_i$ are real-valued signals, what we get are  two circuits each with separate (real-valued) input and separate real-valued output, and internal cross-connections between the two circuits. We have that
\begin{align}
y[n] &= x[n] + a\cdot x[n-1]\\
&\Downarrow\\
y_r[n] + jy_i[n] &= x_r[n] + jx_i[n]
+ \big(a_r +ja_i\big)\big(x_r[n-1]+jx_i[n-1\big])\\
&\Downarrow\\
y_r[n] &= x_r[n] + a_rx_r[n-1]-a_ix_i[n-1]\tag{2}\\
y_i[n] &= x_i[n] + a_rx_i[n-1] + a_ix_r[n-1]\tag{3}
\end{align}
We see that the circuit with input $x_r$ has output $y_r$ that depends on its inputs $x_r[n]$ and $x_r[n-1]$ (as it should) in exactly the same way as $(1)$ but also includes the cross-connection term $-\alpha_ix_i[n-1]$ which depends on $x_i$, the input to the other circuit. Similarly, the circuit with input $x_i$ has output $y_i$ which depends not just on $x_i[n]$ and $x_i[n-1]$ (as it should) in exactly the same way as $(1)$ but also includes the cross-connection term $\alpha_i x_r[n-1]$ which depends on $x_r$, the input to the other circuit. Thus, the block diagram above can be duplicated with some minor changes to labels (e.g. $x[n]$ is replaced by $x_r[n]$ in one circuit and by $x_i[n]$ in the other circuit, and similarly for $x[n-1]$, $a$ and $y[n]$) and the cross-connections added.  I will leave the actual task of sketching the full diagram to those interested in crossing every i and dotting every t.

And that's it.  No messing around with $\sin()$ and $\cos()$ or orthogonality, no phase-shifting by $\pi/2$, or time-shifting


Converting a simple analog Butterworth filter into a digital filter I was just studying an old circuit analysis textbook that was describing how to design a Butterworth filter, and that seemed easy enough.. then, I started to wonder if I can take this analog filter and convert it into a digital filter. Its not really an exercise in the textbook, i was just curious how to convert the analog filter into a digital filter just for fun without the heavy DSP theory.

So I was tried to taking a toy Butterworth filter to do just that.  For example, Let's suppose I had an analog filter:

$$H_a(j\Omega) = \frac{1}{(1+j\Omega)(2+j\Omega)}$$

and i wanted to convert this into a digital filter with say a sampling period $T=200\pi$ rad/sec, and neglecting the effects of aliasing, using this formula:

$$H(e^{j\omega}) = H_a(j\Omega)\Big|_{\Omega = \omega/T}$$

What would $H(z)$ and $h[n]$ look like for the digital filter?

Converting the analog filter $H_a(s)$ into a digital filter $H_d(z)$ using the bilinear transform where T is the sampling period:

$\Large H_d(z) = H_a(s)\bigg|_{s=\frac{2}{T}\frac{z-1}{z+1}}$

Example:

Given a first order Butterworth filter

$H_a(s) = \frac{1}{1+RCs}$

$H_d(z) = H_a\bigg( \frac{2}{T} \frac{z-1}{z+1} \bigg)$

$H_d(z) = \frac{1}{1+RC\Big(\frac{2}{T}\frac{z-1}{z+1} \Big)}$

$H_d(z) = \frac{1}{1+RC(\frac{2}{T}\frac{z-1}{z+1})}$

$H_d(z) = \frac{1+z}{(1-2RC/T)+(1+2RC/T)z}$

$H_d(z) = \frac{1+z^{-1}}{(1+2RC/T)+(1-RC/T)z^{-1}}$

The coefficients of the denominator are the 'feed-backward' coefficients and the coefficients of the numerator are the 'feed-forward'  coefficients used to implement a real-time digital filter.


Normalize speech features per dimensions or across dimensions I am working on pipelines for audio processing (audio features that will go into a neural net to create some embeddings).

We are using two features : Filter banks, or MFCC.

I was wondering what was the good approach with normalization :

Should we normalize those frames on each dimension (each frequency / mfcc coefficient), or normalize across all dimensions ?

This could lead to important differences in the relative importance of features .

Maybe a good starter would be to aggregate those features across your frames into one single feature vector, compute the mean and standard deviation and perform z-score normalization.


FSK: non-integer vs integer rmodulation index, peculiar peaks I've been experimenting with frequency shift keying (FSK) and found a peculiar pattern I cannot explain.
I noticed that there is a significant change to the shape of the power spectral density of the fsk signals, only depending on their modulation index.

For readability I split the discussion into the two cases integer and non-integer modulation index:

Integer modulation index:

The two characteristical FSK flanks are symmetrical and they  show very strong peaks in their middle. The distance between those two peaks is exactly the shift of the FSK signal.  The image below nicely shows those peaks and the symmetrical FSK flanks in the PSD and STFT of one FSK signal.

Non-integer modulation index:

The two characteristical FSK peaks are not symmetrical, there is considerable skewness to be observed.The peaks found in integer modulation index, are not to be found. The image below shows a similar non-integer modulation index example:

I wonder how that difference in the PSD of the signals comes to pass?

With integer modulation indices, each symbol is actually an integer number of full oscillations. You DFT that, and get a sharp, discrete spectrum of tones, convolved with the pulse shape.

With non-integer indices, you get "cutoff" oscillations. That is,  we need to dsitinguish between two cases:

Continuous-Phase FSK: your next oscillation starts with a phase determined by the previous symbol(s) (as the phase continues where you left off)
Non-continuous-Phase FSK: your symbol starts with a specific phase, no matter what the end phase of the last symbol was.

In the non-continuous case, you get a PSK-alike spectral component (because you suddenly change the phase). If you've been (virtually) letting the subcarrier oscillator run through while the other one was active, you'd only need convolve that PSK spectrum with dirac impulses to form your spectrum. If you "stopped" it, and always started a symbol with the same phase, you'd even get the PSK effect if constantly sending the same symbol. If you then also send both possible symbols alternatingly, you'd have to convolve that with the abrupt symbol phase change, so you'd get a slightly biased PSD. I think this is what we're seeing here.

EDIT Found out this is CP-FSK:

Well, the fact that you don't have full oscillations definitely explains why you don't get diracs in spectrum – especially since you're not really observing the PSD (which is an "invisible" property of a stochastic process), but a PSD estimate, very likely based simply on a mag² of a DFT – in other words, oscillations that don't fit into the DFT size an integer amount of times cannot ever be sharp spikes. In reality, this poses no problem (because for reliable detection, you need signal power, not peak height, and the power within the "broadened" peaks is, thanks to Parseval, the same).

Non-symmetry is an interesting effect here, and I wouldn't have expected it to be so clear (it should have been there, by the pure fact that the support of your subcarrier shapes in frequency domain would not be limited to "their" half of all frequencies, just not as clearly). The fact that the "right" half of the spectrum seems to be consitently higher worries me a bit – is it possible there's channel model/analog reasons for this?


physical meaning of fourier transform I have a question regarding the fourier transform. If I plot the frequency against the fourier transform for a periodic signal and I get a peak, What is the physics behind it? I want to know the physical meaning of it.

Why we get a peak and what does it signify?

Thanks

I read your question a little differently than the other people who have answered.  The role of the Fourier transform in physics is as a bi-directional one-to-one mapping between conjugate variables.  The most famous example of a conjugate variable pair is position/momentum, but there are many others. Some examples of conjugate variable pairs:

position and momentum over Planck's constant $[x,p/\hbar]$
time and energy $[t,E]$
aperture and the sine of the diffraction angle divided by the wavelength
$\left[x,\frac{\sin\theta}{\lambda}\right]$

One of the consequences in quantum mechanics is that the operators (quantum mechanics uses differential operators which are related to the classical variables) do not commute.  So, using the position/momentum example; $xp\neq px$.  This leads to an uncertainty relation between conjugate variables.

The usual statement of the Heisenberg uncertianty principle is that the more well one knows the position of a quantum particle, the less one can know about the momentum.  Mathematically;
$$
\sigma_x\sigma_ p\geq\frac{\hbar}{2}.
$$
However, the more general statement of the uncertainty principle is that the uncertainty between any pair of conjugate variables is given by
$$
\sigma_1\sigma_2\geq\frac12.
$$
This uncertainty relation leads to some rather strange effects in quantum physics.  For example, the uncertainty relation between energy and time allow for the violation of energy conservation on very brief timescales, a law which is sacred in classical physics.

If you are looking for more information on the uses of the Fourier transform in physics and engineering with many worked examples, have a look at this free textbook by J. F. James.


What is the sampling rate in this case? A signal $x=\sin(\pi t/4)$
sampled at every $t=1$ sec.
so $T=1$ sec
Sampled signal is then $x=\sin(\pi n/4)$,
What is the sampling rate in this case?,
And according to Nyquist sampling rate what should be the minimum sampling rate to sample this signal?

For the sine wave $\sin ((\pi/4)t)$ , the minimum sampling frequency is greater than 1/4 Hz. That amounts to a sampling period less than 4 seconds. In your case, sampling period of 1 second is therefore a valid choice.


Can we use Tiva C series TM4C123GXL for signal processing experiment purposes? I've newly bought texas instruments Tiva C series TM4C123GXL processing board for developing DSP application project. Is that processing board is good for that purpose? Is that well work with Matlab? Or C2000 Piccolo TMS320F28075 32-bit microcontroller is good for my purpose ?

This question is apparently out of scope. But I would like to provide you a simplest guideline so that you can by yourself decide on whether a given microprocessor would be adequate for various DSP tasks.

1- The Floating Point Capabilities:
Nowadays any processor claiming high quality DSP computation must have strong floating point hardware support. Look for its FLOPs (floating point operations) throughput. Also look for single precision (32 bits) or double precision (64 bits) FP capabilities. For most multimedia purposes (speech, music, video, image) 32 bit is more than enough. For more accurate scientific computations 64 bits is necessary. This being said, the experts can use integer capabilities of a microprocessor to implement most DSP algorithms. But that would require a lower level handling of the mpu. Depending on what kind of an algorithm you want to implement, you should first determine number of floating point operations per second and then check for the MPU if it can support that.

2- Parallel Execution Capabilities (SIMD)
Many DSP algorithms (especially image and video prcessing types, or communication channel processing) will benefit from single instruction multiple data (SIMD) processing pipeline. It can be expensive to include in a simple MPU but if you have this capability, it's worth the money.

3- ADC / DAC capabilities are vital as you will eventually be inputting an analog signal, convert to digital, process it and cpnvert to analog before outputting. So look for for High Speed, High Accuracy ADC and DAC. Or if your application can tolerate low accuracy go with that one.  And there must be enough number of channels so that if you have multiple signals being input or output you can use that.

4- Memory: in addition to limited program code memory for several applications, data has to be buffered inside RAM. this ram can be on-chip or out of chip. First determine how much ram you will need and then look for wheter it can be on chip (very fast) or out of chip.


Why not stop with the first level decomposition in wavelet? Typically, We do wavelet decomposition with the last one level.
But I wonder, why do decomposition withh the last one level?
In the other word, why not stop with the first decomposition level?

The number of levels (frequency range represented) is discretionary to the user usually automatically set as to include a limited number of frequencies which is dependant on the lenght of the series and which is representative of most of the signal.

Applying a wavelet function splits the signal in two coefficients (scaling and detail coefficients) for a number of levels which represent a range of frequencies of this split defined by selecting the number of decomposition levels. Selecting the number of decomposition levels is subjective as per the frequencies you may like to extract.


Circuit simulator with variable simulation time step. How to obtain FFT afterwards? I know that I am not the first one with this problem, but I didn't manage to find the proper answer, so I hope you would like to help me.

From an electric circuit simulator, I simulated my circuit in time domain and saved a signal in a .txt file (which contains both time and signal values).
Based on that .txt file, I would like to produce a FFT to see the signal in the frequency domain.

Unfortunately, the circuit simulator has a variable time step, which gives me a headache to obtain the correct FFT values (two-sided or one-sided is second priority at this point)

I managed to get some FFT out of my data, but I am having difficulties verifying the correctness of it, and hence appreciate any help on this matter. My simulation files give me different values for Average and RMS values, hence the confusion. Thank you.
The link to my .txt file is given below (since I couldn't find a way to upload the txt file here directly):
https://file.io/PXuD9sWR

My Python code:

import numpy as np
from scipy.fftpack import fft
import matplotlib.pyplot as plt
import pandas as pd

buckstruct = pd.read_csv('buck.txt', sep = '\\t', engine='python')
bucktime = buckstruct.iloc[:,0]
bucktime = bucktime.values.tolist()
buckcurrent = buckstruct.iloc[:,1]
buckcurrent = buckcurrent.values.tolist()
bucktime_flt = []
buckcurrent_flt = []

# Filter out the start-up transient. Save only steady state values
for i in range(len(bucktime)):
    if bucktime[i] > 0.002:
        bucktime_flt.append(bucktime[i])
        buckcurrent_flt.append(buckcurrent[i])

plt.plot(bucktime_flt, buckcurrent_flt)
buckfft = fft(buckcurrent_flt)
buckfft_flt = []

# Double the amplitude for harmonics as a first step to converter from two-sided to single-sided FFT
for i in range(len(buckfft)):
    if i == 0:
        buckfft_flt.append(1.0/len(bucktime_flt)*abs(buckfft[i]))
    else:
        buckfft_flt.append(2.0/len(bucktime_flt)*abs(buckfft[i]))

plt.plot(buckfft_flt)

Thank you for your response. I registered to this website, so not sure if my name is the same one or not anymore, but I am the OP.
I'm not quite sure I fully understand the solution on how the resampling and interpolation can help me. If I downsample (let's say every 100th sample), I will have a 100x smaller Dataframe size, but the elements are not evenly separated, i.e. the sampling time between two samples is not constant.

If I upsample afterwards to create empty entries in my Datafra, I get my original Dataframe size. So far, so good.
But the interpolate function somehow seem not to work in my case. I'm new to Pandas, so please be patient :)

import math
import numpy as np
from scipy.fftpack import fft
import matplotlib.pyplot as plt
import pandas as pd

buckstruct = pd.read_csv('buck.txt', sep = '\\t', engine='python')
buckstruct.index = pd.to_datetime(buckstruct.index, unit='s')
buckstruct_downsampled = buckstruct.resample('100s').sum()
buckstruct_upsampled = buckstruct_downsampled.resample('10s').sum()
buckstruct_interpolated = buckstruct_upsampled.interpolate()

EDIT1: Previous contribution merged into this answer:

Thanks for your reply. Yes it's LTSpice, and thanks for the tool. I will look at it. Though it would be nice to also have a solution in Python.

Let me ask the question the other way around, since I am no expert on FFT. Maybe there is no need for this kind of workaround after all.

Given the RAW data samples from my .txt file, i.e. not resampled, interpolated or whatever, and by simply running the FFT, I receive the following FFT result:

So, I have something around 1.75 (y-Axis) at 0 (x-Axis), and something around 0.5 (y-Axis) at 300 (x-Axis).

The 300 is NOT the frequency, because the frequency of my electronic circuit is 100kHz. So I need to find a way to map the 300 in the graph to 100kHz to match my simulations. This is fairly easy if the sample rate was constant. But is it also possible if the sample rate in the .txt is non-uniform?

Thank you,

EDIT2: Solution to the problem: Translating a non-uniformed dataset into a uniformed dataset.

This can be achieved using Scipy's interpolation function. Below is a minimal working example that does not need the initial dataset, but uses 10 consecutive data samples coming from the initial dataset. As one can see, the sample time is not constant which makes it difficult to perform a FFT and to extract the frequencies. Using the interpolation function, constant sample time can be achieved and the resulting data set can be used for FFT calculation. The procedure is to obtain the interpolation function first, and then for a given data length to perform the actual interpolation and to use the interpolated data set for FFT purposes. I hope this minimum working example can be useful for others.
Thank you,

# Minimum working example to show the effectiveness of the interpolation function

import matplotlib.pyplot as plt
from scipy import interpolate
import numpy as np

time_raw = [0.0009999650511607746,
 0.000999965566641282,
 0.0009999660821217893,
 0.0009999672947956213,
 0.0009999682836389868,
 0.0009999696068367931,
 0.0009999715724226046,
 0.000999978044536965,
 0.0009999994033983994,
 0.001]

signal_raw = [1.132772,
 1.132642,
 1.132511,
 1.132205,
 1.1319549999999998,
 1.13162,
 1.131124,
 1.1294879999999998,
 1.124091,
 1.12394]

fxxx = interpolate.interp1d(time_raw, signal_raw)

num = len(time_raw)
xx = np.linspace(time_raw[0], time_raw[-1], num)
yy = fxxx(xx)

plt.figure(1)
plt.plot(time_raw, signal_raw,'bo-', label='Original')
plt.plot(xx,yy,'g.-', label='Interpolated')
plt.ylim([1.120, 1.135])
plt.legend();

The resulting FFT for the initial problem with the actual data set from the simulation is then as expected. My apologies for this messy thread.


mathematical equation for Power Line communication channel How can I derive the mathematical equation for the powerline communication channel to capture the input data signal, the power line channel, and noise?

Agreed to others. You need to narrow down your interest to a frequency band as PLC can be achieved as narrow-band or wide-band. Your channel and noise model would depend on the frequency band you are interested in. For instance, if you are interested in broadband communication (e.g. OFDM), then you need to take transmission/reflection coefficients of the network of interest / topology, attenuation factor (which is a function of distance and frequency) and delay of all the paths in the receive path into account. Of course, the channel is also time-varying and the noise is always assumed to be non-Gaussian.

A pretty common channel model assumed for broadband communication would be as the following:

\begin{equation} \label{eqn_PLN_channel}
H(f)=\sum^L_{i=1}\left\{\prod_{k=1}^O \Gamma_{ik} \prod_{m=1}^P
T_{im} A(f, d_i) e^{-j2\pi f\tau_i}\right\}
\end{equation}
where $\Gamma_{ik}$ and $T_{im}$ denote the reflection coefficient
and transmission coefficient, respectively, $A(f, d_i)$ is the
distance and frequency dependent attenuation and $\tau_i$ represents
the delay of $i$th path out of all the receive paths ($L$) at the point of signal reception.

Derivation of this equation is simple and stems from what's physically happening in the network. Though, note that the channel depends both on frequency and time.

On the other hand, usually there are no closed-form estimates of the noise and the noise is numerically estimated based on a model such as Nakagami-m or Generalized Gaussian Distribution.


Possible to Implement Baseband OFDM on Microcontroller? I am wondering if one can implement OFDM on any available microcontroller or Development board in the present Market. Definitely not looking FPGA solution. Specifically, on any Microcontroller is it possible?

OFDM is just math; an OFDM receiver is just an algorithm. As a microcontroller is a universal computer (i.e. equivalent to a Turing machine), it can execute any algorithm. So, the answer is yes.
Whether or not it is performant enough / has enough memory for the specific OFDM system you have in mind is something nobody can tell you without knowing exactly what kind of OFDM system, what performance bounds and what microcontroller you have in mind.


What is use of a gaussian filter over histogram? I am working in the field of digital image processing and currently i am reading this paper and I have two questions about it....

After applying a gaussian filter to a histogram, will the pixel values of the new histogram be changed or not?
And will the sum of pixels in new histogram be the same or not?

For you questions:
1. After applying gaussian filter on a histogram, the pixel value of new histogram will be changed.
2. The sum of pixels in new histogram is almost impossible to remain unchanged.

Visually speaking, after your applying the gaussian filter (low pass), the histogram shall become more smooth than before. Thus, the new histogram is definitely changed. With regards to the sum of pixels in new histogram, you may consider the problem a little bit differently.

Assume that you have a normalized image histogram $$p(x) = {\# {\rm pixels\,\, of\,\, value\,\,} x\over N},$$where $N$ is the total number of pixels in an image. Then $p(x)$ is also a p.m.f., because $\sum_{x}p(x)=1$. Then the sum of pixels $$S=E[x]\times N,$$
where $E[x]=\sum_{x}xp(x)$ is the mean pixel value over the entire image.

Let's say your have a 1D gaussian filter as
$$f(y)={1\over \sqrt{2\pi}\sigma_y}\exp\left(-{(y-\mu_y)^2\over 2\sigma_y^2}\right).$$ Commonly, people set $\mu_y=0$ as a zero-mean gaussian, and $\sigma_y^2$ is the variance.  As a result, the new histogram after applying the gaussian filter can be written as
$$p'(x) = \sum_zp(z)f(x-z),$$
which can be considered as the gaussian weighted sum of the original neighbor histogram stems. Alternative, you may consider the relationship among the new histogram, the old histogram and the gaussian filter as
$$p'(x) = p(x)*f(x)$$
where $*$ denotes the convolution operator.

This notation then brings out an alternative interpretation from random variables (r.v.). Assume r.v. $a\sim p(a)$, and r.v. $b\sim f(b)$, then r.v. $$c=a+b \,\sim\, p'(c),$$
when $a$ and $b$ are independent. Though it is somewhat hard to believe at the first glance, this interpretation tells you that the theoretical histogram of a noisy image (which is corrupt by the noise following the same gaussian distribution you used in filtering) is the identical to the histogram that you filtering the original histogram with that gaussian filter. One thing that I need to point out here is that even though $c\sim p'(c)$ doesnot mean that you will see the exact same distribution, because all what you can see is realizations of this distribution, but not itself directly.  Please see the following figures for details. As you can see, red and green curves match very well, but they are still slightly different.

Finally, we are ready to answer your second questions whether the sum of pixels changes before and after filtering. Using our previous derivations, you will see this question is equivalent to ask whether
$$E[c] =E[a] ?$$
where $c\sim p'(c)$ and $a\sim p(a)$. If $b$ follows a zero-mean gaussian, then $E[c]=E[a]+E[b]=E[a]$. Otherwise, they are not equal. Once again, you may not expect to see these values are identical, because $E[c]$ is only a theoretical value, and the actual observed mean value can be deviate from this theoretical value (the degree of how much the observed mean deviate from the theoretical one depends on the number of pixels, the more pixels the closer these two means).

Here is some Matlab code that generates the above figures. Play this demo, and you will see what I mean.

I = imread('cameraman.tif');
[p,val] = imhist(I);
p = p./numel(I); % <== p(x) original histogram
figure,plot(p)
sigma = 5;
f = normpdf(-20:20,0,sigma); % <== f(x) gaussian distribution
p1 = conv(p,f); % <== p'(x) new histogram after applying gaussian filtering
hold on, plot(val(1)-20:val(end)+20,p1,'r--','Marker','o')

% alternative interpreation
N = randn(256)*sigma; % <== white gaussian noise following p.d.f f(x)
I2 = N+double(I); % <== a noisy observation by adding noise to the original image
[p2,val] = hist(I2(:),val(1)-20:val(end)+20); % <== p''(x) the hisogram of the noisy image
p2 = p2./numel(I);
hold on, plot(val,p2,'g-.','Marker','*')
legend('original histogram', 'histogram after applying gaussian filter', 'histogram after adding gaussian noise')

figure,
subplot(131),imagesc(I),axis image,colormap(gray),title('original image I');
subplot(132),imagesc(N),axis image,colormap(gray),title('white gausian noise N');
subplot(133),imagesc(I2),axis image,colormap(gray),title('noisy image I2 = I+N');


Inverse wavelet transform on wavelet coefficients I have a data set which consists of Morlet coefficients for 7 channels in 18 bands. When I try to use inverse DWT in Matlab to get the raw data in time, using this code:

X = idwt(Data(:(#band),:(#time),1(#channel)),[],'morl');

I get this error :

***********************************************
ERROR ...
-----------------------------------------------
 wfilters ---> The wavelet morl is not valid!
***********************************************

any Idea what should I do ? I'm pretty noob so any explanation is great help for me

Because, as far as I know, there is no discrete wavelet scheme related to the morlet wavelet, admitting neither orthogonal nor biorthogonal implementation. See for instance the discussion in Complex Morlet function and DWT.

Using waveinfo('morl'), you indeed get that: DWT: no

 Information on Morlet wavelet.
    Morlet Wavelet
    Definition:
    morl(x) = exp(-x^2/2) * cos(5x)
    Family                  Morlet
    Short name              morl

    Orthogonal              no
    Biorthogonal            no
    Compact support         no
    DWT                     no
    CWT                     possible

    Support width           infinite
    Effective support       [-4 4]
    Symmetry                yes

So you might obtain something only (and I am not sure it is valid for your data) using some form of inverse continuous wavelet transform, for instance in Matlab  icwt.


Accelerometer - coordinate system transformation I'm getting some accelerometer readings from an Android phone, but it comes in on the phone's coordinate system. I want to apply a transformation to put the acceleration in the world coordinate system (e.g. so Y will always be up and down, regardless of how the phone is rotated)

Android allows you to get the accelerometer values. It can also calculate a rotation matrix (using a gyroscope and a magnetic sensor), so all I need to do is multiple the acceleration vector through by this matrix and I have transformed values. I followed the steps here: https://stackoverflow.com/questions/14963190/calculate-acceleration-in-reference-to-true-north

However, there are 3 sensors at play here: acceleration, gyro, magnetic. The way android operates is event based, so whenever it detects a change in any of these sensors, it will give you back a value.

Given these are 3 different sensors, the events for each one is very likely to come in at different times, and you may get more readings from one sensor than another within some time period

My question is, does varying amounts of spacing between readings for the 3 different sensors matter for doing coordinate system transformations? And if so, what is the best solution for doing this type of transformation with event based sensors?

It will call your function each time any motion event gets fired. This can affect the timing of css animations responding to sensor input.  Try getting the orientation values and updating your transformations from within an interval.


CIE-Lab changes on Light intensity changes I am planning to use CIE-Lab color model for an image segmentaion method.

If the light intensity of sun, some led source changes on the objectm then does it affects the a and b value of the data also?

And if different color of light source is focused in the object, that affect has it on the lab values of given image?

I hope I get it right. You want to reduce effect of the light sources and perform the image segmentation. So basically you want to perform color based image segmentation.
If you do want to do that on a variety of data the first step should be performing the white-balancing technique on your images, and adjust the color temperature on the same values. This ensures the same type of your illumination source, even if the sources change over time (e.g. color temperature of daylight changes). The type of the light source definitely affects your color values. Unfortunately this can be performed only on RAW image datatypes, but if you expect change of light sources and you are keen on accurate results, I would strongly recommend the white-balance process.

For image color segmentation I do recommend mentioned CIELAB color space, because of it's independence. It does not depend on the input devices as well. The big advantage here is the possibility of separation the Lightness channel and perform segmentation only on the color ones (a*, b*). That is why this color space is more useful than for example RGB, which includes information of lightness in it's R, G, B components and the results are significantly biased due to that. However the L channel of the CIELAB corresponds more to the psychometric lightness rather than the real luminance.  The problem might be the intensity of light source, EV as well. If you get saturated and also really dark parts of the image you will not be able to restore the right color values. Other settings of the camera might have impact on the colour values as well. So if you want seriously accurate measurements and results from image color segmentation it might be really tricky. On the other hand, with a "normal" conditions and removing L channel in CIELAB color space, you can reduce the effects of the light sources on colour values as much as possible and there are not many other ways to do that. So ideally take your images with the same settings, white-balance them, transform them to CIELAB, remove the L channel and than perform color based segmentation.


Parallel Goertzel filters in C# I have 2 frequencies in my application: $12000\mathrm{Hz}$, $13000\mathrm{Hz}$.
I want to set a parallel Goertzel filters and every 100ms I want to take the bigger magnitude for detecting the tone received.

Is there any way to implement this?  That is, how can I evaluate two functions in parallel in C#?

I work on WPF, C#

Thanks

First of all, a link to an example of the Goertzel algorithm.

Second, the goertzel algorithm is so fast that you don't need to get into multithreading.

Just create a class that does the goertzel algorithm for a frequency and sampling rate that you define when you create an instance.
Make a list or array of instances of that class, then for every sample, run through that array and have it calculate the goertzel for that sample.

In your case, you just need two instances so you don't really need an array of instances.
Fast and simple.
I've done this very thing with several hundred instances on a fairly old PC, and updating the GUI causes more of a load than the processing.


Is there an audio processing library/API that contains all standard features of a digital audio workstation? I need an audio processing library/API (in any programming language, preferably Java or Python) that is able to do things to sound waves such as: sum/mix sources together, apply EQ, dynamic EQ where you can change a frequency peak or amplitude setting over time, pitch shift, distortion, and convolution reverb, all programmatically without any need for a user interface.

Someone asked "If there were, why would there be a market for premium commercial DAW software?" Obviously this wouldn't be for the actual music-making composer or sound engineer because it would be unusable (it wouldn't have a GUI). I need it to automatically manipulate sound so I can churn out millions of random not-necessarily-musical sound files, without a GUI.

This is not same as the "OpenCV for audio" question, because I don't need intelligent audio feature extraction; I just need the standard features of DAW such as EQ and convolutional reverb, available as an API call, like soundFile.applyEq(myEq).applyConvReverb(myReverb);

For python, I'd try out:

pydub for basic audio manipulation and arrangement (slicing, concatenation, fading, repeating, etc).
For effects processing, pyo is a good python signal processing module with a musical focus. Here's a comprehensive review, with inline clips.

Together you should be able to get good coverage of most DAW functions. There are also fantastic python audio libraries like librosa, in case you do decide to do any feature extraction or audio analysis. Good luck!


Is the phase difference between recorded sounds from left and right ear the same as the phase difference between their HRTF? Pardon my ignorance if my question is silly. Lets say I did a binaural recording and found the phase difference between right and left channel. If the recording was done in anechoic chamber, where there is no reverberation, can I say that this phase difference is same as the phase difference between their HRTFs which they were filtered by (on their way from speaker to my ear)?

Consider a signal with a particular frequency at phase $\alpha$. Passing this through the two filter paths changes the phase of the frequency to $\alpha + L$ and $\alpha + R,$ where $L$ and $R$ are the values of the phase frequency responses of the two HRTF filters at the particular frequency. Yes, the phase difference $(\alpha + R) - (\alpha + L)$ is equal to $R - L,$ which you call more simply the "phase difference" between the HRTF's.


How to use values from the output of a block in GNURadio I want to use the output values of the top_block as a way of directing my flowdiagram. So I want the flowdiagram to take one path if the output of my processing block is 0 and take a different path if it is non-zero. The if self.blocks_multiply_xx_0.out_sig != 0: statement is the condition I am having trouble with.
I have implemented an if statement within the top block which disconnects and connect parts of the diagram accordingly.

I want the if statement to look at the output of the self.blocks_multiply_xx_0 block and send it to a file sink if it is non-zero, or otherwise send it to a null sink.
I want this to be done for each and every value coming out of the self.blocks_multiply_xx_0 block.

However to do so I need a way of accessing that blocks output value.

How can I go about accessing that value, and assessing which path to take based upon that value, each and every time the self.blocks_multiply_xx_0 block produces an output value?

Is the location of the if statement correct?

Will it check each and every processed sample value or does it have to be moved to a different method?

If so where could it go?
from gnuradio import gr
from gnuradio import audio
from gnuradio import analog
from PyQt4 import Qt
from gnuradio import blocks
import threading
import time
class my_top_block(gr.top_block):
def init(self):
gr.top_block.init(self)
     samp_rate = 32000
     amp = 0.1
     # Blocks
     ##################################################
     self.blocks_multiply_xx_0 = blocks.multiply_vff(1)
     self.audio_sink_0 = audio.sink(samp_rate, '', True)
     self.blocks_file_sink_0 = blocks.file_sink(gr.sizeof_float*1, 'C:\SDR\GNU files\speed_test_SSD\speed_test_SSD.iq', False)
     self.blocks_null_sink_0 = blocks.null_sink(gr.sizeof_float * 1)
     self.analog_sig_source_x_0_0 = analog.sig_source_f(samp_rate, analog.GR_SIN_WAVE, 350, amp)
     self.analog_sig_source_x_0 = analog.sig_source_f(samp_rate, analog.GR_SQR_WAVE, 0.5, 1)

     ##################################################
     # Connections
     ##################################################
     #connect square wave with multiply
     self.connect((self.analog_sig_source_x_0, 0), (self.blocks_multiply_xx_0, 0))
     # connect sine wave with multiply
     self.connect((self.analog_sig_source_x_0_0, 0), (self.blocks_multiply_xx_0, 1))
     #connect multiply output with audio sink
     self.connect((self.blocks_multiply_xx_0, 0), (self.audio_sink_0, 0))
     #connect multiply out to with null sink
     self.connect((self.blocks_multiply_xx_0, 0), (self.blocks_file_sink_0, 0))
     #connect multiply out to with null sink
     self.connect((self.blocks_multiply_xx_0, 0), (self.blocks_null_sink_0, 0))

     #if its non-zero send to file sink, otherwise send to null sink
       if  self.blocks_multiply_xx_0.out_sig != 0:
           self.disconnect(self.connect((self.blocks_multiply_xx_0, 0), (self.blocks_null_sink_0, 0)))
           self.connect((self.blocks_mgetinultiply_xx_0, 0), (self.blocks_file_sink_0, 0))
       else:
         self.disconnect(self.connect((self.blocks_multiply_xx_0, 0), (self.blocks_file_sink_0, 0)))
         self.connect((self.blocks_multiply_xx_0, 0), (self.blocks_null_sink_0 ,0))

def d():
my_top_block().run()
def n():
t0 = time.time()
time.sleep(5)
t1 = time.time()
t = t1 - t0
print 'time = ', t
def ram2disk():
#take the contents in written to Q: drive
#write these to a new file in C: drive
#execute this in a new thread occurring simultaneously as the write to ram

if name == 'main':
t = threading.Thread(name='non-daemon', target=n)
d = threading.Thread(name='daemon', target=d)
d.setDaemon(True)

d.start()
t.start()</pre>

I want to use the output values of the top_block

A top_block has no output; a signal processing block in a flowgraph has.

So I want the flowdiagram to take one path if the output of my processing block is 0 and take a different path if it is non-zero.

Ah, you kind of want to "demux" something. So write a simple demux block that has two inputs and multiple outputs; copy from the first input to any of the outputs, depending on your second input.

I have implemented an if statement within the top block which disconnects and connect parts of the diagram accordingly.

Don't do that. Flow graph reconfiguration is asynchronous to signal processing in GNU Radio, i.e. this won't work, and especially not on a per-sample basis.


Order of applying two Butterworth filters (distributivity) I have a function $f(a, b)$ that applies a linear transformation on a signal $a$, given another signal $b$.

I would like to filter the result of $f$ using a 2nd order Butterworth low-pass filter ($\mbox{butter}()$).

My question is, would first applying the filter on the individual signals give the same result as applying the transformation on the raw data and filtering afterwards?

In other words, does the following hold: $\mbox{butter}(f(a, b)) = f(\mbox{butter}(a), \mbox{butter}(b))$ ?

PS. would this property be called distributivity? I can't really figure this out..

The short answer: it depends.

The longer version: From MIT OpenCourseWare notes:

  Specifically, if we have several LTI systems cascaded together, the
  output generated by an input to the overall cascade combination does not depend
  on the order in which the systems are cascaded.

In your specific case, you have a cascade of two blocks: one which is a 2nd order Butterworth low-pass filter (which can be shown to be a linear-time-invariant (LTI) system), and another which is your linear transformation function $f(a,b)$. The question then boils down to whether $f(a,b)$, in addition to being a linear transformation, also happens to be a time-invariant transformation.

If $f(a,b)$ happens to be an LTI transformation, then the cascade $\mbox{butter(f(a,b))}$ is equivalent to $f(\mbox{butter}(a),\mbox{butter}(b))$. Thus, in that case one could say that $\mbox{butter}$ is distributive over $f(a,b)$.

If $f(a,b)$ isn't time invariant, then we can't make that claim as to whether the result depends on the order of operations.


Random demodulator in matrix form I went through a paper in random demodulator. It states that If the sampling rate is $R$ ($R$ is less then $W$, $W=$band limit of a signal in Hz), and assume that $W$ is divided by $R$. Then each sample is the sum of $W/R$ consecutive entries of the demodulated signal. Does it mean that the signal is compressed by taking $W/R$ samples at a time?

The random demodulator (in this definition) compresses the signal by taking samples that are each linear combinations of W/R Nyquist rate samples. So it does not take W/R samples at a time; it takes one sample at a time, each of which contain information corresponding to several (W/R) Nyquist rate samples. It is then up to a reconstruction algorithm to subsequently "disentangle" these into the original signal (approximately) at Nyquist rate.

The random demodulator works by "smearing" the content of the input signal all over the spectrum (by multiplication with the spreading sequence). Provided that the input signal is sufficiently sparse, we can now reconstruct it from just a small portion of the resulting spectrum (selected by the low-pass filter).


What does this signal filter do? (coming from finance) In finance, many "indicators" are used in order to make prediction about if a price has an increasing or decreasing trend over time.
Among them, MACD is a popular one. There are many slighlty different definitions (standard rolling mean vs. exponential rolling mean, etc.), but one of them is :

the price is $x[n]$   (let's say the sampling rate is 1 second, but we could also use 1 day...)
$y_1[n]$ is the 12-period rolling mean, $y_2[n]$ is the 26-period rolling mean, i.e. :

$$y_1[n] = \frac{1}{12} x[n] + \frac{1}{12} x[n-1] + ... + \frac{1}{12} x[n-11]$$ and the same for $y_2$.

the MACD is defined by $y_1 - y2$, i.e. :

$$macd[n] = \frac{7}{156} x[n] + \frac{7}{156} x[n-1] + ... + \frac{7}{156} x[n-11] - \frac{1}{26} x[n-12] - \frac{1}{26} x[n-13] + ... - \frac{1}{26} x[n-25]$$

Question:

How to have a signal processing interpretation of what's going on with $macd[n]$? I can see it's a linear filter, but how to know what it is doing? Can we draw a frequency response?
It also seems that, because of rolling mean, some lag / delay is introduced.

Remark: When people use MACD, they don't stop with $macd[n]$: they define $signal[n]$ as the 9-period rolling mean of $macd[n]$:

$$signal[n] = \frac{1}{9} macd[n] + \frac{1}{9} macd[n-1] + ... + \frac{1}{9} macd[n-8]$$

and the indicator which is finally used is  (this sounds like a magic trick! why these choices?):

$$h[n] = macd[n] - signal[n] =
-\frac{14}{351}x[n]-\frac{49}{1404} x[n-1]-\frac{7}{234} x[n-2]-\frac{35}{1404} x[n-3]-\frac{7}{351} x[n-4]-\frac{7}{468} x[n-5]-\frac{7}{702} x[n-6]-\frac{7}{1404} x[n-7]+\frac{2}{27} x[n-12]+\frac{7}{108} x[n-13]+\frac{1}{18} x[n-14]+\frac{5}{108} x[n-15]+\frac{1}{27}x[n-16]+\frac{1}{36}x[n-17]+\frac{1}{54}x[n-18]+\frac{1}{108}x[n-19]-\frac{4}{117}x[n-26]-\frac{7}{234}x[n-27]-\frac{1}{39}x[n-28]-\frac{5}{234}x[n-29]-\frac{2}{117}x[n-30]-\frac{1}{78}x[n-31]-\frac{1}{117}x[n-32]-\frac{1}{234}x[n-33]$$

When $h[n]$ is positive the trend of $x[n]$ is more or less increasing, and when $h[n]$ is negative, the trend of $x[n]$ is more or less decreasing.

Question:

What does $h[n]$ do, with a DSP point of view?

I thought you stock market tech analysts always used moving averages (rolling means) or "exponential moving averagers" for your smoothing operations. Anyway, your $macd[n]$ process is neither of those. As far as I can tell the frequency magnitude response of your $macd[n]$ process looks like the following:

Your $macd[n]$ process does not smooth its input data in the normal sense. It amplifies input data that has a periodicity (repetition rate) of 33.3 days per cycle. In the above plot you can clearly see the "mainlobe" that Fat32 was referring to. If your closing price data contains a periodicity of roughly 33 days then the output of your $macd[n]$ process will be a large value. Otherwise, the output of your $macd[n]$ process will be a small value. I don't know if such a frequency response makes stock market technical analysis sense, or not.

As for your $h[n]$ process, I believe its frequency magnitude response looks like the following:

Your $h[n]$ process appears to be a process that amplifies input data that has a periodicity (repetition rate) of 21.9 days per cycle. In DSP we use processes like your $macd[n]$ and $h[n]$ for filtering when we want to detect the presence of components in our input signal that have periodicities of 33.3 or 21.9 time units per repetition. We call such processes "resonators" or "bandpass filters."


From Orthogonality to DTFT I have met this question, but cannot prove it using DTFT definition.
Given: $g$ is a discrete sequence filter and:
$$ g \in l^2(Z)$$
$$\langle g_n, g_{n-2k} \rangle = \delta_{k}$$
Prove:
$$|G(e^{j \omega})|^2 + |G(e^{j(\omega + \pi)})|^2 = 2$$

Could you guys give me some direction on this one?
Thank a lot

Here a few hints to get you started:

Define a sequence $x[k]=\langle g_n, g_{n-k} \rangle$ and write it as a convolution. From that it should be easy to find its DTFT $X(e^{j\omega})$ in terms of $G(e^{j\omega})$.
Next note that the given sequence is just a downsampled version of $x[k]$, so its spectrum is an aliased version of $X(e^{j\omega})$. That's how you end up with those two terms on the left-hand side of the frequency-domain representation.


Where can I get passport size face photo database for my image classification experiment? I am just trying to classify images using kNN algorithm with the aid of MATLAB. So now I am in need of some collection of images of peoples' faces. Where can I get some open database of faces (both females and males).

Here is an up-to-date index of face image databases commonly used by researchers.

http://www.face-rec.org/databases/

It contains links to the common datasets such as The Yale Face Database (regular and "B"), The Color FERET Database, The AT&T Database of Faces, etc.

It should cover your needs.


How does reconstruction function relate to watershed segmentation? I'm trying to understand how watershed work by reading The Morphological Approach to Segmentation: The Watershed Transformation S. Beucher and F. Meyer

This document mentions the reconstruction function. But I don't understand how reconstruction function relates to the watershed algorithm?

More detail, in reconstruction algorithm implementation, it is implemented by ordered queue (below image)

And input of algorithm :

First, I understand parcel and film in function but I have no idea about it in image. Can you give me some example about it on image and how can we get film known parcel ?

Finally, result of this algorithm :

I know that is not result of watershed, so how does it relate to watershed's output ?

The reconstruction function corresponds to the "flooding" of a basin.

The relationship between the reconstruction and the Watershed Transform is as follows:

Given a minimum seed point, that is, the minimum grayscale value on some image, the watershed algorithm starts "flooding" the image. Intuitively, "flooding" is an operation of assigning a label to all the pixels that are adjacent to the seed point and have a higher grayscale value than the seed point.

Notice these two points here "adjacent" (to what? How to assess adjacency?) and "higher grayscale value" (Progressive "flooding" is essentially tracking the circumference of a basin at some height from its lowest point).

On the first step of this process, we start with a seed point. A seed point has 4 or 8 neighbours and in this case it would be very easy to decide which of those neighbouring pixels should be "flooded" in the next step.

But, as this process evolves, what started as a "seed point" now becomes a surface whose outer boundary can meander in any possible way. This "surface" is the surface of the flooded basin.

Mathematically, this is expressed on page 446 by:

$$W_{i+1}(f) = \left[IZ_{Z_{i+1}(f)}(X_i(f))\right] \cup m_{i+1}(f) $$

Where:

$W$ is the surface of the watershed
$IZ$ is the "Zone of Influence" of a section
$Z$ is the lower section of the image $f$
$f$ is the image itself (it is viewed here as a two dimensional function that maps a two dimensional set of integers to a grayscale value, that is $f : \mathbb{Z}^2 \rightarrow \mathbb{Z}$ or perhaps the more familiar $f(m,n)$ that would return the grayscale value of $f$ at point $m,n$).
$X$ is the upper section of the image $f$
$m$ is the set of local minima of $f$ at some level $i$.

The "Zone Of Influence" of a single point or more generally, of a region of points of interest is the set of points at a finite minimum geodesic distance from the region of interest. This is specified on page 438, in two steps, as:

$$z_X(Y_i) = \left\{ x \in X : \begin{matrix} d_X(x,Y_i) \text{ finite} \\ \forall j \neq i, d_X(x,Y_j) < d_X(x,Y_j) \end{matrix} \right\}$$

$$IZ_X(Y) = \bigcup_i z_X(Y_i)$$

Where:

$z$ is the Zone of Influence
$d$ is the geodesic distance between two points
$x$ is a point in the region of interest
$X$ is the region of interest
$Y$ is the rest of the points in a specific image

BUT, as you can probably observe from the discussion so far all these quantities of interest are binary or come as the result of operations between binary objects. Remember that a grayscale image is some $f: \mathbb{Z}^2 \rightarrow \mathbb{Z}$, so, how do we calculate that $d$ which is then taking part in these "selection" operators?

Distances are calculated between sets of points (including point-to-point, expressed as a set of one element). A grayscale image contains numbers within the range defined by the pixel depth, how do we establish "distance" between any two of them?

The "trick" (sic) to converting a grayscale image to a binary one is via the "sections". The sections are the first thing that is defined on page 434 as upper section $X_i(f) = \left\{ x \in \mathbb{Z}^2:f(x) \geq i \right\}$ and lower section $Z_i(f) = \left\{ x \in \mathbb{Z}^2:f(x) \leq i \right\}$.

This is a view of an image as being composed of slices (sections) and looking roughly like this:

With this view in mind, it is now possible to compute geodesic distances between sections because they simply represent sets of points and therefore also possible to compute the reconstruction of a section. The reconstruction of a section refers to all the pixels that are at some finite distance from the section. This is the $R_X(Y) = \left\{ x \in X: \exists y \in Y, d_X(x,y) \text{finite} \right\}$ that is defined on page 437 and to a large part it drives the definition of the Zone Of Influence $IZ$.

The crucial point to note here is that if you derive $R$'s repetitively for each section $i$ of the image, these $R$'s will be "nested" and "tangential" to each other allowing them to be used to simulate this effect of "flooding" from lower levels to higher levels. This is explained on page 439 from "Consider two functions..." to "...The dual operation is a closing.".

In conclusion, given a grayscale image $f:\mathbb{Z}^2 \rightarrow \mathbb{Z}$, we define the upper and lower sections at grayscale level $i$ as $X_i$ and $Z_i$ respectively. Given a section and a metric of distance $d$, we define the reconstruction function $R$. The reconstruction function $R$ helps in defining the Zone of Influence $IZ$ which is what the Watershed Transform "tracks", over a number of iterations to perform the segmentation.

Hope this helps

(Image retrieved from google images)


Matlab: What is the proper way to calculate mean square error for complex numbers I am confused regarding the calculation of mean square error involving complex numbers. Considering, the true channel coefficients to be

h_true =

   0.7071 + 0.7071i
   0.4243 + 0.4243i
   0.2121 + 0.2121i

and the estimates to be

 h_e1 = 0.8100+0.8100i
 h_e2 = 0.5100 + 0.5100i
 h_e3 =  0.1200 +  0.1200i;

For training epoch 1, I calculated the error between the estimates and the actual  as

e = ((h_e1 - h_true(1)).^2) + ((h_e2 - h_true(2)).^2) + ((h_e3 - h_true(3)).^2)   ;

Then the mean square error is

mse = e/3;

which gives zero for the real part and some small number for the imaginary part.

mse = 0 + 0.0176i

Should I be calculating the error separetely for iamginary and real part or is this okay?

The mean square error is

e = ((abs(h_e1 - h_true(1))^2) + (abs(h_e2 - h_true(2))^2) + (abs(h_e3 - h_true(3))^2))/3;

But it is tediuos!

What if you had $1000$ terms instead of $3$? Do you want to add one by one?

Do it in the vector form for more clarity and ease of implementation:

$$\mathrm{mse}=\mathsf{E}(|\mathbf{h}-\hat{\mathbf{h}}|^2)$$

In MATLAB, you can use the function mse()

 h_true = [ 0.7071 + 0.7071i
            0.4243 + 0.4243i
            0.2121 + 0.2121i ];

 h_e = [ 0.8100 + 0.8100i
         0.5100 + 0.5100i
         0.1200 + 0.1200i ];

 e = mean(abs(h_e-h_true).^2); % or eqivalently e = mse(abs(h_e-h_true));


Sudden loudness changes producing unwanted noises? What's that effect called? Is that one of the reasons for having attack and release parameters when implementing synths? And is there a generic method to avoid such sound artifacts?

Example: I have implemented a simple synth that uses a midi control channel to adjust loudness by applying the 128 discrete MIDI values directly to the generated audio signal. Whenever one touches the MIDI controller wheel, pretty nasty, scratch-like sounds get produced.

You are changing the volume in discrete steps, applied at random points in the audio signal.

If the signal happens to be close to zero, you won't hear much of anything.

If the signal is close to a peak, then you will get a "step" change in the signal.  This audible.

When you roll the volume up and down that way, you get a kind of "zipper" sounding noise.

What you want to do is to apply volume changes only when the signal is close to zero.  The "zippers" still occur, but are much less annoying - and may be completely inaudible if you change volume only at times when the signal is at zero.

Looking at the code you referred to in a comment, it appears you aren't changing volume when the amplitude is close to zero.

You are changing volume when the signal is less than zero.  What you want to do is to change volume when the absolute value of the signal is less than some small value.

Start with like 0.1, then try 0.01.  You should notice far less "zipper noise" that way.

Sort of like this:

if (abs(amplitude) < 0.1d) {
   currentVolume = targetVolume;
}


Recommended resources for noise reduction I know the question What Resources Are Recommended for an Introduction to Signal Processing (DSP)? and I have already read some general DSP books, such as Rick Lyons's "Understanding DSP,"  and also have browsed through Oppenheim's.

More specifically, what (printed) book would you recommend for denoising techniques in general?

Including techniques adapted for audio, spectral subtraction,  Wiener filtering,  etc. and if possible adapted to discrete signal processing with an algorithm/engineering point of view
  (i.e. rather x[k] than continuous x(t) and integrals...)

I just ordered a used version of Advanced digital signal processing and noise reduction, Vaseghi, but I'm looking for other recommendations.

Signal denosing is a well-studied technique in signal processing. It first began using simple techniques such as filtering. In this approach, the emphasis is laid on designing filters which can perform denosing techniques in a fast and efficient manner.

Later, when wavelet theory was developed, some researchers used wavelets to denosing the signal. Essentially, we can apply thresholding techniques (soft or hard threshold) on the coefficients of the wavelets. On reconstruction using the dominant coefficients, one can denoise the signal.

In the early 2000s, sparsity-based techniques were introduced where the input signal is represented using an over complete dictionary. In the reconstruction step, sparsity inducing norms were used to extract only the dominant columns of the dictionary. Popular approaches include LASSO and total variation denosing.

Reference

Useful reference with MATLAB codes http://www.numerical-tours.com/matlab/


How to represent the nonlinear model as a state space in Unscented Kalman Filter There is an Autoregressive model of order 1 (AR(1)) that is excited by a non-linear signal as the input:
$$x_t = \rho x_{t-1} + u_t \tag{1}$$
The time series $u_t$ is generated from a nonlinear map,
$$u_t = f(u_{t-1},\mathbf{w}) \tag{2}$$ where $f$ is the nonlinear function.
The observations are $$y_t = x_t + v_t \tag{3}$$ where $v_t$ is the measurement noise that is an Additive White Gaussian Noise.

Q1: Can I re-write the model (1) as a state space in the following way:

$$x_t = Ax_{t-1} + f(u_{t-1},\mathbf{w}) $$

$$y_t =  x_t+ v_t \tag{4}$$

Is the above representation correct? If not then I shall be grateful for the correct technique to represent it.

Rather than write

$$
x_t = Ax_{t-1} + f(u_{t-1},\mathbf{w})
$$

as the state update equation, I'd write:

$$
\xi = \left[ \begin{array}{c}
x_t\\
u_t
\end{array}
\right]
$$

and then
$$
\xi_t = g(\xi_{t-1},A,\mathbf{w})
$$

so that
$$
y_t = \left[ 1\ \ \ 0 \right] \xi_t + v_t
$$

Then you could apply the non-additive noise formulation of the EKF to get your filter equations.


What is the best way to represent audio visually? (x-post from UX) Original Question: https://ux.stackexchange.com/q/23040/16006

I've only taken some basic signal analysis courses, so I might be missing some things.

Purely theoretical question:

What methods exist for representing audio?

What methods could be made for representing audio, more specifically musical audio?

So far, I'm aware of:

Viewing the waveform (Soundcloud does this), mostly useless except for seeing "loudness"
Spectral analysis (Example), good for seeing frequency and "loudness"

Essentially I'm wondering if there is a way one could "see" the notes, beats, and so on of a song, visually.

Right off the top of my head I can think of displaying 3 differently colored waves over time representing treble, mid, bass in a soundcloud-like container with the section playing (or moused-over) being magnified, with the surrounding waveforms being compressed into the corners (like a wide-angle lens effect).

EDIT: I don't know where this could be used, this was just born out of my frustration with current audio visualizing technology.

I imagine having a 3d graph of a spectral analysis over time (Ninja Edit: apparently known as spectrogram) would be the "best" solution since you see everything but it might not be the most elegant and it might not be portable to places like soundcloud.

Even current spectrum analysis is hard to decipher (Too low level for images):

FL Studio wave editor

I'm essentially wondering what might work for casual users, and for people wondering ahead of time how the song will play out.

What a human (or their ear-brain) perceives in sound is a psychoacoustic phenomena, and may or may not be exactly related to the actual audio as recorded.  e.g. the exact notes, beats and instruments that a human "hears" may be influenced by visual cues, memory of other similar music, and the musical context around the actual sound of the note in question.


What is Causing a Noise Signal in USB Audio Recording? I am experiencing the following peculiar noise signal when recording from my Behringer Xenyx Control 2 over USB.

There is frequency at around 2293 hz as show:

And over time, it produces the following plot:

What on earth is doing this?

It's NOT USB problem. If you there was anything wrong with the digital connection, you would get drop outs or nothing at all. The error detection/correction on the USB bus is binary: either you got all the bits or you didn't.

This is clearly analog noise, probably from switching power supply. It happens most likely in the analog front end of the mixer. Once it's digitized, it can't pick up this type of noise anymore.

You can try a different USB recording device or throw a USB hub in the path to isolate the iMAC power from the Xenyx power see if it helps. It's also entirely possible that the Xenyx will do this no matter what: Behringer has a reputation for great value, but also for sloppy design and poor quality control. The Xenyx is a blatant ripoff of the Mackie Onyx (as the name implies).


How to check my CRC of my RF packet? The "problem" is that I want to validate if my received RF packet is valid or not valid, in other words I want to validate if my RF packet corrupted or not corrupted .
the transmitter that Im using is micro controller launchpad called c1350 , so through the datasheets for the CC13xx I discovered ( in page 5, section 4) that it uses CRC32 (page 704, section 8.2.3.8). So, the implementation of crc that my transmitter cc13xx useis CRC32.... my question is:
I've the recieved packet and I know its CRC in terms of binary bits, so now how can I check by CRC if my recived packet is valid or not? my packet before transmitting "before getting to the air" wasn't having CRC .. and once transmitted "on the air" it has CRC (received packet has CRC).
what approach/algorithm should I do/implement in order to check if my received packet corrupted or not?
CRC implementation is CRC32 .
by datasheet cc13xx: enter link description here we discover at (page 5, section 4)
It uses CRC32: enter link description here (page 704, section 8.2.3.8).

There doesn't seem to be much discussion of CRC in this forum but on the sister site math.SE, there are several postings. See, for example, CRC Computation or What is the theory behind seeding a CRC.


Fast Hartley Transform Implementation in MATLAB I want to implement Fast Hartley Transform (Specifically Discrete Hartley Transform) in a script file in MATLAB. Does anyone know have a reference implementation of this in MATLAB or another language so that a complete novice can understand it?

You can resort to the properties of the DHT with the DFT, given in the wikipedia article. For a real-valued function $f$, its DHT in terms of the DFT is given by

$$
\mathcal{H}f = \Re\{\mathcal{F}f\}-\Im\{\mathcal{F}f\}.
$$

So, a function for fast Hartley transform can be given by

function Hf = fht(f)
F = fft(f);
Hf = real(F) - imag(F);
end

Though, you need to take care of the normalization factors, and tdhis is up to you and how you define the scale factors of the DHT.


SDR : Signal divided between the edges of the receiving window I'm using a SDR (specifically the RTL-SDR) to analyse the spectrum and i noticed a phenomena that does not seem logical.

The signal I want to capture is centred on 868.3 MHz and the central frequency of the SDR is centred on 869.3 MHz with a bandwidth of 2 MSamples. What I observed is that the 2 halves of the signal are displayed on both sides of the receiving window instead of just the half between 868.3MHz and 868.3MHz - 0.125MHZ ( as my signal bandwidth equals 125 KHz ) .

What do we call this phenomena and what is the explanation behind it ?

The first image is when the central frequency of the SDR is 868.3MHz

The second image is when the central frequency ofthe SDR is 869.3MHz

That's aliasing in its most classic form.

The anti-alias filter in your receiver isn't perfect, so that frequencies slightly outside of what you want to sample appear on the other end.


What is Local Mean Filter? The research paper "Multidirectional Scratch Detection and Restoration in Digitized Old Images" says that,

  4.1. Preprocessing. The preprocessing step aims to enhance image features along a set of chosen directions. First, image is grey-scaled
  and filtered with a sharpening filter (we subtract from the image its
  local-mean filtered version), thus eliminating the DC component.

Now, I am Googling the term "Local Mean Filter" but there is nothing available like that.

Can anyone please provide me any reference of "Local-Mean Filter"?

What filter are they using for Sharpening?

They probably just wanted to say that the image was blurred (by some local method, e.g. convolution with a Gaussian kernel) in a more scientific way.

On the sharpening: They don't sharpen the image directly. What they do is to blur the image and the subtract the blurred version from the original. The result is the same as sharpening directly via an appropriate filter.

Why this works

Assume the image $I$ and an blurring kernel $g$ (which could be a Gaussian kernel). In Fourier-Space, blurring the image by convolution becomes a multiplication:

$\mathscr{F}(g\star I) = \mathscr{F}(g)\mathscr{F}(I)$

with $\mathscr{F}(\cdot)$ being the Fourier-Transformation.

Then $I - g \star I$ becomes: $\mathscr{F}(I) - \mathscr{F}(g)\mathscr{F}(I) = \mathscr{F}(I)(1-\mathscr{F}(g))$.

Since $\mathscr{F}(g)$ enhances the low-frequency components in in Fourier-Space, $1-\mathscr{F}(g)$ enhances the higher frequencies (I assume here that $\mathscr{F}(g)$ is normalized appropriately in magnitude, which can be achieved via scaling factors). So, in image-space, the resulting image is a sharpened version of $I$.


Concept About Estimated Standard Deviation I am looking for the concept about how to estimate standard deviation.
Actually I'm not sure how can I get a concept the estimate standard deviation ?
If you know the concept, then would you let me know it here?
Also if you have any reference or examples, it will be help to me.
p.s. Is this different between "estimate standard deviation" and "estimated standard deviation" ? I think it seems to the same.
Updated
What does different between "Method of Least Squares" and "estimate standard deviation"?
Are those the same concept? I don't know estimate standard deviation.

Given data $ { \left\{ {x}_{i} \right\} }_{i = 1}^{N} $ the Empirical STD of the data is well defined:

$$ STD = \sqrt{ \frac{1}{N - 1} \sum_{i = 1}^{N} { \left( {x}_{i} - \bar{x} \right) }^{2} } $$

Where $ \bar{x} $ is the empirical mean of the data given by:

$$ \bar{x} = \frac{1}{N} \sum_{i = 1}^{N} {x}_{i} $$

Now, if there's a model on the data (Such as signal and noise with certain CDF) the empirical calculation should be altered accordingly.

For instance, given a signal which is linear with AWGN the STD of the noise can be estimated by removing the linear estimated signal first.


How to differentiate the product signal $f(t)\theta(t)$, where $\theta(t)$ is Heaviside's unit step function? What is the derivative (in the engineer's sense) of the causal function $f(t)\theta(t)$, where $\theta$ is the Heaviside unit step function?

I've seen the formula $f'(t)\theta(t)+f(0)\delta(t)$, where $\delta$ is Dirac's delta. This looks like a kind of "product rule": differentiating the product gives $f'\theta+f\theta'$, but $\theta'$ is $\delta$, and $\color{blue}{f(t)\delta(t)=f(0)\delta(t)}$.

If this is right, I don't understand the following argument, from the solutions manual to Oppenheim and Wilsky's Signals and Systems. The solutions manual says the derivative of the function $2e^{-3t}\theta(t-1)$ is

$$-6e^{-3t}\theta(t-1)+\color{red}{2}\delta(t-1)$$

It's the second term I don't understand. Using the "product rule" heuristic, the second term should be $2e^{-3t}\delta(t-1)$, which using the blue formula above gives $\color{red}{2e^{-3}}$ times the delayed delta function, not just twice the delayed delta function.

Is the solutions manual wrong?

(Cross-posted on MSE.)

The answer has already been given in Fat32's comment: if the solution manual says what you claim then it is wrong and you are right.

What I want to add here is that the product rule for distributions is by no means some "engineering heuristic" or some dubious magic, but it can be proved in a rigorous way (see e.g. this document, p. 5) by treating the combined function $f(t)\theta(t)$ as a distribution. I point this out because you don't seem to be sure about this (you mention: kind of "product rule", "product rule" heuristic).

Furthermore, also the fact that the (generalized) derivative of the step function equals the Dirac delta impulse can be shown in a rigorous way. So even if we use notation that seems appalling to most mathematicians, we engineers tread on pretty safe ground when it comes to the use of distributions as given in your example.

A very good introduction to distributions and singularity functions that is rigorous yet accessible for engineers can be found in The Fourier Integral and Its Applications by A. Papoulis.


Confusion about PSD and RMS Let's say I have a noise power-spectral-density (PSD) which is not flat and ranges from 0 to $f_1$ Hz in frequency. As we know, the total area under the PSD is equal to the total average power of the signal (according to Parseval's theorem). Let's say this value is $\sigma^2$. This means the value of the auto-correlation function at zero-time (which is the total power in the signal) will also be equal to $\sigma^2$.

Now assume we want to approximate the above colored noise with a band-limited (from 0 to $f_1$) white Gaussian noise with the same noise power (i.e. $\sigma^2$). That means the variance of the AWGN is $\sigma^2$, the auto-correlation function would be a delta function (i.e. $\sigma^2 * \delta(t)$) or more precisely a sinc() function, and the PSD would be flat (equal to $\sigma^2$) from 0 to $f_1$. The total area under the PSD should be equal to the total average power of the noise. The total area under the PSD in this case simply would be $\sigma^2 * f_1$ which is not equal to the initial assumption of the total power $\sigma^2$. What am I missing here? Should I consider the flat PSD value is $\sigma^2/ f_1$ instead of $\sigma^2$ (however this does not satisfy the auto-correlation vs PSD relation)?

I think what your missing is that the Power Spectral Density, as a density is the power per unit frequency, so assuming units of Hz as is typically done, $\sigma^2/f_1$ in your case is the power per Hz; and therefore to get the total power, $\sigma^2$, you would multiply that by the total frequency $f_1$ resulting is the expected $\sigma^2$.


What is wrong with RMS and Maximum delay spread Matlab code? pls, help check this matlab code. The code is to determine the RMS and Maximum delay spread. After running the code below, I was not getting desirable results, the Max. delay spread is giving the same value in all the data except where i have more than two spikes. I also got some negative max. delay spread in some cases, I guess I did do one or two things right. Pls, can someone help to check the work-ability of this Matlab code below? Thanks in anticipation.
this is the matlab code:
clc

clear

format short

ts= 1e-06;  % ts=1

T = readtable('Rx8.csv');

kk=table2array(T);

z=length(kk);

allData=zeros(z,2);

num_5=0;

sum_p_i=0;

Mean_delay=0;

Mean_t2bar=0;

for ii=2:z
% if(ii)

s(ii,1)=kk(ii,1)+1i*kk(ii,2);

pdp=abs(s(ii,1))^2; % Values for PDP

allData(ii,1)=pdp; % Store PDP in matrix
%

tau=(ii-1)*ts; % delay_n

allData(ii,2) = tau; % Store delay_n in matrix
%
num_5=num_5+(pdp*tau); % Equation (2.18)

sum_p_i=sum_p_i+pdp
%
t2bar=(pdp*tau^2)/pdp; %nth_delay_square_bar

end

tbar=sqrt(num_5/ sum_p_i);

for jj=2:z

tau_rms= sqrt(( allData(jj,1)* (allData(jj,2)-tbar)^2)/sum_p_i);
end

[aaa, bbb]=max(allData(:,2))

[aaa, bbb]=max(allData(:,1))

[B,I] = sort(allData(:,1),1,'descend')

delay1=allData(I(2),2)

delay2=allData(I(1),2)

% delay2=allData(I(3),2)

maxdelay=delay2-delay1

tbar

tau_rms

plot(allData(:,2), allData(:,1))

% RMS_delay=sqrt(Mean_t2bar-Mean_delay^2)

%
% t_pdp=1:length(pdp)

% plot(s(:,2),pdp)
%
% [a, b]=max(pdp)

% [a2, b2]=min(pdp)
% plot

% Amp1=sum(abs(s.^2))

% Amp2= rms(s)^2

% Amp3=norm(s)

Not an answer but some guidance how to post this type of questions which is too long for a comment (@mods: feel free to delete if that's inappropriate).

We typically don't debug code. Focus on algorithmic questions. If possible ask about equations, not the implementation of the equations in code.
Define your terms (like "maximum delay spread") using an external reference or preferably with an equation to make sure we are all talking about the same thing.
If you use an equation (which you should), define all symbols used in the equation.
If you post code, make sure it can be run. We don't have your data file so we can't really reproduce your problem.
Do basic debugging first yourself: Run your code with the most simple signal with a known answer. Once that runs ok, add more complicated signals ( but still with a known answer). Only after this is all working, proceed to unknown signals and answers
If you still have a problem, please post why you think this is. What is the answer you expect and why? What is the answer that you are getting and how is different from your expectation? Why does the problem doesn't seem to appear on the basic debugging signals you have used in step 5.