Why exactly is a bit/radix reverse required when calculating the fft? When implementing the fft, you always include a bit reverse at the end of the process to put the data in the correct order. It seems to be because the fft proceeds by decimating the input array into odd/even sections (for radix-2 algorithms, for other algorithms e.g. prime factor you have to do a radix reverse), the final stage of the algorithm involves multiplications (by several exponential factors) on the individual array members. These can be labelled by a series of EOEEO etc (E = even, O=odd) depending on the process of decimation and whether the subarray was an odd/even decimate. Why does replacing E with 1 and O with 0 result in the bins requiring a final permutation. This has always bothered me about the fft.

One reason typical implementations of a radix-2 FFT require a permutation at the beginning or end is to reduce the use of any temporary memory when using a procedural/sequential processing engine.  When using a single complex array, larger than, say, the CPUs register set, for both the input, output and intermediate results, the intermediate butterfly results have to go somewhere in that array.  But the place most appropriate for a non-permuted result might still be occupied by unprocessed input data.  So you could either use more memory to save the unprocessed inputs somewhere else, use more memory to put the intermediate results somewhere else, or just put the results in an array location where the input data is no longer needed (e.g. the butterfly input locations).

If you allocate a full new array for each FFT radix-2 layer, no post-shuffling will be required.  You could also ping-pong 2 arrays, but that's still double the memory.

If you had enough parallel ALU units that no temporary result are required, and the entire FFT result is available at once at the ALU outputs, then you could also write out the FFT result without any permutation.


Detection milliseconds signal using a DVB-T dongle as SDR in C (need idea) I have probably amateurish problem. Assume that we have a USB tuner DVB-T + connected antenna.

I would like to write a program that will detect the signals of at most a few milliseconds.

Signal looks classically as follows:

Someone can help me how to do it? I want to write it in C.

Edit: I do not need, of course, the final code. I need ideas.

Edit-2:

I add details of the signal:

Frequency: approx. $850 MHz$ (may be more)

The measurement of dispersion: $> 250 pc.cm^{-3}$

The length of the signal: $> 0.01 ms$

Peak Flux: (I'm not sure what value we can get)

Your C program can link to librtlsdr (source available on github), set the frequency and sample rate, then filter the incoming IQ sample stream from the RTL USB dongle (demodulate if necessary), and (one possibile algorithm out of many) implement an envelope matched filter.


Minimum period of a signal What is the minimum period P (in samples) of the signal $e^{j(\frac{M}{N} )*2πn}$ for the following values of M and N?

M=1,N=3
M=5,N=7
M=35,N=15

I have got the answer for the first pair of values, which is 3. But I could not get the other two correct. Please help.

It's already answered in the comments so just to close it out: the answers are 3, 7 and 3. The periodicity is basically given by the denominator of the fraction in front of 2*pi*n after the fraction has been simplified (no common divisors between nominator and denominator).  The simplified fractions for the cases are 1/3, 5/7 and 7/3 so the periodicity is 3, 7 and 3 respectively.


Estimation of step amplitude using the CUSUM algorithm I used the CUSUM algorithm to detect steps in data. Basically the data looks like this, the data has a constant amplitude and then there is a rapid variation or a step. For example, the signal has a value of 1000 and rapidly drops to 950 or 900. With the CUSUM algorithm, I can detect easily these changes and the time at which they occur. But now I would also like to automate the amplitude estimations of these steps. Eventually, I'll need to implement this in real-time ( a small delay is tolerable for the amplitude estimation)

One solution I would like to try is this one :

take N samples before the step and take N samples after the step and subtract the respective averages.

Is there a better solution?

Edit :
- The transition from 1000 to 900, for example, might take 10 samples in all.

Your straightforward idea of using moving averages might actually work quite well. Since you ask "is there a better solution", I must ask what you mean by "better."

Change point detection is a well studied problem in statistics and signal processing. The problem has many variations depending on your optimality criteria, and prior knowledge about the data. For example: do you care about catching the change point as quickly as possible? Or do you care about the precise location of the change? Or do you care about the exact amplitudes before and after the change, or just the difference in amplitudes?

Since you have a real-time requirement, I would search for "online Bayesian change point detection" which will show you lots of papers, and many packages in R to play with.

https://www.r-bloggers.com/a-simple-intro-to-bayesian-change-point-analysis/

https://www.jstatsoft.org/v23/i03/paper

https://hips.seas.harvard.edu/files/adams-changepoint-tr-2007.pdf

Detection of Abrupt Changes: Theory and Application, Basseville and Nikiforov


Select talk not music from podcast I listen to a podcast with interesting talk separated with horrible music.
I'd like to rewrite the podcast automatically to delete this "noise".

Edit: the talk and music is not simultaneously, it's talk for a couple of minutes, then a horrible song, then talk etc.
Erring on the side of caution is OK, some seconds of the beginning and end of music is OK.

Ideas?

If this is a typical case of a podcast that sounds like a radio broadcast then you could discriminate between music and speech purely based on the "shape" of the spectrum.

I am not sure if there is special software that could handle this but I will provide a simple example using Audacity further below.

The main idea is that a mixed track coming from a variety of different instruments would be expected to occupy a wider range of the spectrum than just someone talking.

Here is a quick test using Audacity:

Just by visual inspection, a set of "noisy bands" appear in the spectrum view. These are the locations where it is either music played along with someone speaking or just music playing.

My settings for this image (Edit->Preferences->Spectrogram): Window-size 4096, window:Blackman, minFreq:0, maxFreq:8kHz, Gain:20, Range:80, freqGain:0, use grayscale.

To switch to this view, load your podcast into Audacity and then click on its track name and select "Spectrogram".

To achieve the same thing in an automatic way, you would have to create a very simple classifier. The classifier would break down the incoming signal into "chunks" and obtain their spectrum. This is roughly what Audacity does to generate each one of the vertical columns that make up the spectrogram of the sound. The classifier would then "scan" these columns from left to right deciding what each column represents ("Music", "Speech", ..., other?). The decision would be based on the relative amplitude of each harmonic in the spectrum. After this, you could decide to drop or silence those frames that have been characterised as "Music".

But all this means writing some form of code in a programming language (e.g. Python). I can expand this response with the simplest of classifiers that could work for this problem but only if that would be seen as useful (?).

Of course, this technique would work with mixed tracks whose full spectrum (including all their possible harmonics) extends well beyond the spectrum of the human voice. There might indeed be some special cases as mentioned earlier by robert bristow-johnson where the discrimination is much more difficult and more complex techniques would be required but for a generic podcast, just basing the decision on the spectrum might do just fine.

By the way, the podcast I used in this example is this one.


Is Hilbert transform not defined for complex signals? Is Hilbert transform not defined for complex signals?  In MATLAB, the function hilbert ignores if you give a complex sequence as input.  Why?

The Hilbert transform can be applied to complex functions of a real variable. E.g., the Hilbert transform of the complex exponential $e^{j\omega_0t}$, $\omega_0>0$, is given by

$$\mathcal{H}\{e^{j\omega_0t}\}=-je^{j\omega_0t},\qquad\omega_0>0$$

The problem you encounter has to do with Matlab's implementation of the function hilbert.m. It is designed for real-valued input sequences and it will ignore any imaginary part. Note that despite its name this function does not simply return the Hilbert transform of the input vector, but it computes the corresponding analytic signal, i.e. it returns a complex vector, the real part of which is equal to the input vector, and the imaginary part of which is the Hilbert transform of the input vector.

So if for whatever reason you want to compute the Hilbert transform of a complex vector, you need to do the following:

x = ...    % some complex vector
xr = real(x);
xi = imag(x);
xr_ = imag(hilbert(xr));    % Hilbert transform of real part
xi_ = imag(hilbert(xi));    % Hilbert transform of imaginary part
x_ = xr_ + 1i * xi_;        % Hilbert transform of complex vector x


What is the importance of signal bandwidth in radar systems for velocity and localization estimation So I know about the basics of doppler shift which relates the velocity of a moving object with the frequency carrier shift seen at the reciever.
But, I don't understand the importance of the transmitted signal bandwidth to the process?
I can use the same the same carrier frequency but with different BW...
Any clarification please?
Thanks

Bandwidth affects the "gate/cell width" of your measurements. If you think about an LFM waveform reflecting back off an object and going through a matched filter to determine the range, a signal with a very low bandwidth is going to have a more ambiguous return time (range) than a more aggressive chirp. The rough calculation for the smallest range resolution for a certain bandwidth is $\frac{c}{2B}$.
I imagine the relationship is reversed for doppler processing. A signal with low bandwidth would provide a lot of doppler resolution because there isn't a great deal of ambiguity as to the return's center frequency, so the doppler gate width should increase proportional to bandwidth. That is a guess, though.
Radar and other signal processing disciplines involve these tradeoffs a lot. Certainty in frequency versus certainty in time and vice versa, maximum  unambiguous range versus improvement factor from pulse repetition frequency, etc.


Examples of Independent and uncorrelated data in real-life, and ways to measure/detect them We always hear about this vector of data VS this other vector of data being independent from each other, or uncorrelated, etc, and while it is easy to come across the math regarding those two concepts, I want to tie them into examples from real-life, and also find ways to measure this relationship.

From this stand point, I am looking for examples of two signals that are of the following combinations: (I will start with some):

Two signals that are independent AND (necessarily) uncorrelated:

The noise from a car engine (call it $v_1[n]$) and your voice ($v_2[n]$) as you are talking.
A recording of humidity every day ($v_1[n]$) and the dow-jones index ($v_2[n]$).

Q1) How would you measure/prove that they are independent with those two vectors in hand? We know that independence means that the product of their pdfs equals their joint pdf, and thats great, but with those two vectors in hand, how does one prove their independence?

Two signals that are NOT independent, but still uncorrelated:

Q2) I cant think of any examples here ... what would some examples be? I know we can measure correlatedness by taking the cross-correlation of two such vectors, but how would we prove that they are also NOT independent?

Two signals that are correlated:

A vector measuring an opera singer's voice in the main hall, $v_1[n]$, while someone records her voice from somewhere inside the building, say in the rehearsal room ($v_2[n]$).
If you continuously measured your heart rate in your car, ($v_1[n]$), and also measured the intensity of blue lights impinging on your rear windshield ($v_2[n]$)... I am guessing those would be very correlated... :-)

Q3) Related to q2, but in the case of measuring cross-correlation from this empirical stand point, is it enough to look at the dot product of those vectors (since that is the value at the peak of their cross-correlation)? Why would we care about other values in the cross-corr function?

Thanks again, the more examples given the better for building of the intuition!

A few elements... (I know that this is not exhaustive, a more complete answer should probably mention moments)

Q1

To check whether two distributions are independent, you need to measure how similar their joint distribution $p(x,y)$ is to the product of their marginal distribution $p(x) \times p(y)$. To this purpose, you can use any distance between distributions. If you use the Kullback-Leibler divergence to compare those distributions, you will consider the quantity:

$\int_x \int_y p(x, y) \log \frac{p(x, y)}{p(x) p(y)} dx dy$

And you will have recognized... the Mutual Information! The lower it is, the more independent the variables are.

More practically, to compute this quantity from your observations, you can either estimate the densities $p(x)$, $p(y)$, $p(x, y)$ from your data using a Kernel density estimator and do a numerical integration on a fine grid ; or just quantify your data into $N$ bins and use the expression of the Mutual Information for discrete distributions.

Q2

From the Wikipedia page on statistical independence and correlation:

At the exception of the last example, these 2D distributions $p(x, y)$ have uncorrelated (diagonal covariance matrix), but not independent, marginal distributions $p(x)$ and $p(y)$.

Q3

There are indeed situations in which you might look at all the values of the cross-correlation functions. They arise, for example, in audio signal processing. Consider two microphones capturing the same source, but distant from a few meters. The cross-correlation of the two signals will have a strong-peak at the lag corresponding to the distance between microphones divided by the speed of sound. If you just look at the cross-correlation at lag 0, you won't see that one signal is a time-shifted version of the other one!


Is it feasible to try and classify an enclosed space as being a. Very crowded b. Less crowded & c. Sparsely crowded using image processing? The idea is to get three classifications for crowd: High, low and medium based on an image captured inside a train compartment. Something like this:

This gets classified in the 'High' category

http://www.spiceflair.com/wp-content/uploads/2012/05/14a.jpg
This perhaps in the 'Medium' & so on.

This sounds feasible (and fun!), but it will be hard to get very accurate.  My guess is that, like your tags suggest, this would be a combination of edge detection / space detection to get a measure of density and face/hand detection to get a measure of human content.  You could then train up an SVM or RandomForest model on the data, tweak, and repeat.


Understanding Convolution Take:

  $$
(u*v)(k) = \sum_{i=-\infty}^\infty u(i)v(k-i).
$$

  The $k$ is there, it's because you want to define
  $$
\ldots\ldots, (u*v)(-3), (u*v)(-2), (u*v)(-1), (u*v)(0), (u*v)(1), (u*v)(2), (u*v)(3), \ldots\ldots
>$$
  etc.  The number in the parentheses is $k$.  Thus, for example, when $k=4$, we have
  \begin{align}
(u*v)(4) = \sum_{i=-\infty}^\infty u(i)v(4-i)
\end{align}
  $$
= \cdots\cdots+u(-3)v(7)+u(-2)v(6)+u(-1)v(5)+u(0)v(4)+u(1)v(3)+u(2)v(2)
>$$
  $$
\phantom{={}} {}+u(3)v(1)+ u(4)v(0) + u(5) v(-1) + u(6)v(-2)+u(7)v(-3)+u(8)v(-4)+\cdots>\cdots.
>$$

Could you give examples of what the 'k' value could stand for?
Can it only be time offset?

Could you give an intuitive explaination why you need to multiply u(-3) [a point left from the y-axes] with v(7) [a point right from the y-axes]! And not u(-3) with v(-3) for example? With a diagram for example.

I think the reason you are getting confused is that you are only considering the convolution of two signals, whereas convolution is not usually between signals but between a signal and an impulse response.

The best way to visualise it is probably with a filter.  For example, in

we can see that the input comes in, and is multiplied by $b_0$ to give its contribution to the output $y$.  Then in the next time step, that input value is multiplied by $b_1$ to give its contribution to the output, and so forth.  In this case your $k$ variable specifies which time in the output sequence you want.  The output signal is given by $y[n] = (x * b)[n]$.

For example, consider the case of $n=4$ (that's the same as $k=4$ in your example).  That means we're looking at the 4th time step (since time = 0).  So the value at the input is $x[4]$, the fourth input value; the value to the right of the first $z^{-1}$ block is the previous input value $x[3]$; the value to the right of the next $z^{-1}$ block is $x[2]$; and so forth.  The output $y[4]$ will be given by
$$y[4] = x[4]b_0 + x[3]b_1 + x[2]b_2 + x[1]b_3 + \cdots + x[-3]b_7 + \cdots$$
 which you should be able to see is the discrete convolution sum, and includes the product of $x[-3]$ with $b[7]$.

Any linear filter or transfer function can be implemented as a convolution, where the output is given by the convolution of the input with the system impulse response.  Try reading up on convolution in a text book, and look particularly at the sliding-strip method of calculating it, which might help with your visualisation.


Interpolation in Contrast Limited Adaptive Histogram Equalization I have been trying to implement the CLAHE algorithm and came across this page which states step by step procedure for the algorithm.

I understand the initial steps to perform HE of tiles in the image. But can't figure out how to perform the interpolation in the last step to remove the artificial boundaries between the tiles? Can someone explain this to me? It would be really helpful if someone explains it with a figure.

In step 8 of the reference you provide, you obtain 4 possible values for a given pixel (based on the values obtained from its 4 neighbours, which I assume are the usual North-East-South_West neighbours).
Then, to avoid introducing tile artifacts, you can use any sufficiently regular interpolation method to obtain the value in the center of these 4 pixels.

For simplicity and computational issues, I suggest that you try bilinear interpolation. It will not introduce ringing or ripple artifacts and will not overflow the dynamic range of the input pixels.


Spectrogram of a chirp and its FFT This question follows from an answer to this other question: Conceptual question on FFT and chirp signal

I wrote a code starting from the spectrogram to compare it with the FFT result.

% Chirp generation
t = [0 0.5 1.0 1.5 2.0]; % time breakpoints
f = [0 200 100 150 300];  % instantaneous frequency breakpoints
p = polyfit(t,f,4);
dt = 0.001;
t = 0:dt:2;
y = chirp(t,p);
L = length(y);

% FFT
NFFT = 2^nextpow2(L);
Y = fft(y,NFFT)/L;
Fs = 1/dt;
f_plot = Fs/2*linspace(0,1,NFFT/2+1);

% Spectrogram
figure()
specgram(y,128,Fs,128,120); % 2D time-freqency display

[S,F,T,P]=spectrogram(y,128,120,128,1E3);
figure()
surf(T,F,10*log10(P),'edgecolor','none')
axis tight

This is what I obtain.

My qyestion is: shouldn't the spectrogram and the FFT have the same amplitude variation? Whereas the spectrogram does not seem to have amplitude variation (the red line is red at the same way everywhere).

Second question, why in the FFT I have high peaks at $100$ and $200 Hz$ only?

FFT calculates spectrum over the whole duration of your signal (integrates energy at each frequency over the time of accumulation) while spectrogram shows frequency variation in time (instantaneous frequency). Divide the whole frequency band of interest in some equal intervals and evaluate how much time signal's frequency falls in each interval. These time periods aren't equal. So by FFT frequencies corresponding to larger time periods gains more power and will be major and vice versa. Now lets look at your spectrogram. We see that it's about 100 and 200 Hz when the most stationarity in frequency exists, there are maximums of the spectrum. Frequencies between 100 and 200 Hz appears 2 times, so we see raising of the spectrum at this span. Low and high frequencies appears only 1 time and frequency at these regions varies faster: the spectrum after FFT decreases.


How to calculate the upsampling rate (sample per symbol) for a recorded raw IQ .wav file? I am not familiar with raw IQ wave files. For tetra-1 signal, the bit rate is $36$ k bps (or $18$ k as the symbol rate). If the Fs saved in the .wav file is $48000$ Hz, what should be the SPS (sample per symbol)? Remember that if SPS is 4 for QPSK, then both I&Q components need to be upsampled by $4$. In this case, should it be calculated as $(48000/2)/18000$?

It is a unit conversion. You have the sampling rate given, $F_s$, and it has units of $\frac{\text{samples}}{\text{second}}$. You also have the symbol rate, $F_{sym}$, and it has units of $\frac{\text{symbols}}{\text{second}}$. You are trying to figure out what is the $\frac{\text{samples}}{\text{symbol}}$ so it is this division:
\begin{align}
\frac{\text{samples}}{\text{symbol}} &= \frac{ F_s }{ F_{sym} } \\
&= \frac{ \big(\frac{\text{samples}}{\text{second}} \big)}{ \big(\frac{\text{symbols}}{\text{second}}
\big) } \\
&= \frac{\text{samples}}{\text{symbol}} \checkmark
\end{align}


Difference between sticky channel Poisson-repeat channel I don't quite understand the difference between a sticky channel and Poisson-repeat channel. A sticky channel allows duplications of symbols sent
over the channel, then the block structure at the sender and receiver will be the same (Taken from "A Survey of Results for Deletion Channels and Related
Synchronization Channels" by Professor Michael Metzenmacher). But I don't understand what is meant by block structure. Is it just referring to the structure the packet, such as the overhead, the payload, etc? A Poisson-repeat channel is defined as: the sender sends a sequence of n bits, and the channel independently replaces each bit with a number of copies (or repeats) that has a discrete Poisson distribution with mean 1. Taken from this blog. But it seems as if the two channels are the same. What exactly is the distinction?

The Poisson-repeat channel is a binary channel where each transmitted bit is repeated at the channel output anywhere from $0$ to $\infty$ times. The repetition distribution is Poisson.

A "sticky channel" is a repeat channel where deletions are forbidden; that is, every transmitted bit is repeated from $1$ to $\infty$ times at the channel output, according to some specified distribution.


Digital filter coefficients from low-pass to high-pass Given I have coefficients a0, a1, a2, b1, and b2, defining the difference equation for a digital filter as:
y[n] = a0 * x[n] + a1 * x[n - 1] + a2 * x[n - 2] - b1 * y[n - 1] - b2 * y[n - 2]
Which defines a low-pass filter with particular cutoff frequency, how can I obtain the coefficients A0, A1, A2, B1, B2, which similarly define a high-pass filter with the same cutoff frequency? I'm aware there are so-called "bandform transformations" for converting a prototype low-pass into a high-pass, but to my knowledge, these are not directly applicable to discrete-time/digital filters, so I am unaware of any way to apply them to this problem.
If these coefficients are derived from a complex-conjugate pair of zeros and/or of poles given the discrete transfer function for the low-pass filter would be:
H(z) = (z - Zero[0]) * (Z - Zero[1])/[(Z - Pole[0]) * (Z - Pole[1])]
Is there then a way to transform this function to the corresponding high-pass filter I'm looking for to get the poles and zeros from the new transfer function?

You can apply a so-called all-pass transformation to a discrete-time low-pass prototype filter in order to convert it to other standard filters (such as high-pass, band-pass, and band-stop). This is accomplished by transforming the complex variable $z$ in the transfer function of the prototype filter by a function $G(z)$ which satisfies $|G(e^{j\omega})|=1$, i.e., $G(z)$ is an all-pass function. This makes sure that the transformation maps the unit circle onto itself, i.e., the frequency response of the new filter is just a shifted and/or warped version of the prototype frequency response.
The most straightforward way to transform a discrete-time low-pass filter to a high-pass filter is to use the trivial transform $G(z)=-z$, i.e.,
$$H_{HP}(z)=H_{LP}(-z)\tag{1}$$
If $H_{LP}(z)$ is given by
$$H_{LP}(z)=\frac{\displaystyle\sum_{k=0}^Nb[k]z^{-k}}{\displaystyle 1+\sum_{k=1}^{N}a[k]z^{-k}}\tag{2}$$
then $H_{HP}(z)$ becomes
$$H_{HP}(z)=\frac{\displaystyle\sum_{k=0}^N(-1)^kb[k]z^{-k}}{\displaystyle 1+\sum_{k=1}^{N}(-1)^ka[k]z^{-k}}\tag{3}$$
This transformation shifts the low-pass frequency response by $\pi$ (i.e., by half the sampling frequency). Consequently, if $\omega_c$ is the cut-off frequency of the low-pass prototype filter, the cut-off frequency of the resulting high-pass filter is given by $\omega'_c=\pi-\omega_c$.
Other cut-off frequencies can be obtained by applying the more general lowpass-to-highpass transformation
$$G(z)=-\frac{z+\alpha}{1+\alpha z},\qquad |\alpha|<1\tag{4}$$
The simple transformation shown above is obtained from $(4)$ by the choice $\alpha=0$.
This and other frequency transformations applicable to discrete-time filters are treated in some detail in Chapter $7.4$ of the third edition of Oppenheim and Schafer's Discrete-Time Signal Processing.
Also take a look at these related questions and their answers: Q1, Q2.


Find the relation between two noisy periodic signals of equal length [Thank you, @leftaroundabout for your helpful comments to my related question. I'm trying to ask a less confusing question here.]

$0\leq\phi\leq 2\pi$ is a parameter that theoretically describes a system. $\alpha=\sin(\phi)$ and $\alpha' = \cos(\phi)$ are observables.

In a measurement, $\phi(t)$ is a continuous periodic function of the time $t$. $\alpha$ can be directly measured and $A_i$ is a noisy over-sampled digital representation of such a measurement. $\beta$ is a different observable, the relation between $\beta$ and $\phi$ shall be studied. $B_i$ is a sampled synchronously with $A_i$ and noisy.

Three questions are interesting and the answer shall be found from $A_i$ and $B_i$:

What is $\delta = \beta(\phi=-\pi/2) - \beta(\phi=+\pi/2)$? This the difference of the $\beta$ values corresponding to the extrema of $\alpha$.
Is $\beta = a * \alpha + c$ a valid model and what is the coefficient $a$?
Can we give a confidence interval for $a$ or $\delta$?

Inspection of the data shows that a better model for $B_i$ is $\beta' = a * \alpha + a' * \alpha' + c$. Does this affect the recipe to find $\delta$ or $a$?

A third alternative model is $\beta'' = f(\phi) + c$. How to reliably find $\delta$?

The nature of the noise is unknown. The noise spectrum has both frequencies less then the modulation frequency of $\phi$ and higher frequencies. A Fourier filter (e.g. lowpass) renders $A_i$ and $B_i$ reasonably smooth.

All right, we can cope with this one.

So first, we should like to have some knowledge about the relation of  $i$ and $\phi$,
$$
  \phi\colon\ \mathbb{Z} \to S^1,\qquad i\ \mapsto\ \phi_0+\frac1{\nu_{\!s}\cdot\omega_0}i
  =: \varphi + \lambda i\quad {}_{^{(\operatorname{mod}2\pi)}}
$$
Where we can essentially forget about the $\operatorname{mod}2\pi$, just view it as $S^1=\mathbb{R}/\!{\sim}$.
There are multiple ways to determine the most likely values of $\varphi$ and $\lambda$; what I'd do is a least-square fit in those variables and $A_M$, perhaps also an offset $A_O$ such that
$$
  \zeta := \sum_i \bigl(A_{\!M}\cdot\sin\,\phi(i) + A_O - A_i\bigr)^2
$$
is minimized. A fitting algorithm like this can be somewhat slowish if you want it to work reliable, I don't know if that's a problem.

Once we have this relation $\phi(i)$, we can just read out the $B_i$ for which $i$ is close to $\tfrac\pi2$ resp. close to $-\tfrac\pi2$. We can simply average over all next-neighbours to these real values, this way we're quite immune from nonlinear response, but rather more subject to noise than with using some kind of interpolation. You want to do it this way if you have a very large amount of data, spanning many periods of $\phi$. If you don't, better use some interpolation between the $B_i$ with $\phi(i)\approx\pm\tfrac\pi2$. I'd try a Friedrichs function as convolution kernel.


Estimate Delay of a Known Signal Delayed by Sub Sample Resolution Given a known signal $ x \left( t \right) $ and its delayed version $ y \left(t, \tau \right) = x \left( t - \tau \right) $.
Both are sampled by Sampling Frequency $ {F}_{s} $ to generate the signals $ x \left[ n \right] $ and $ y \left[ n \right] $.

Assume they are sampled on a time support much larger than their non zero time (Namely, the whole part of the signal where it is not vanished is sampled).

Given those 2 sampled signals, using Cross Correlation (Matched Filter) the Delay Parameter $ \tau $ can be estimated.

Yet, its estimation is limited to the grid defined by the sampling Frequency which means if $ \tau < {T}_{s} = \frac{1}{{F}_{s}} $ the estimation isn't accurate (Actually, the part which is the modulo of the time interval can't be estimated).

Is there an effective method to estimate the delay in Sub Sampling resolution?

Thank You.

Since the Matched Filter is used the Cross Correlation becomes "Auto Correlation" which is assured to be "Symmetric" relative to its maximum.

Hence a good approximation would be to approximate the area around the sampled peak by a Symmetric Function as a parametric model of the function.

The most trivial and easiest model for such function would be the Parabolic Function which is illustrated in the following.

Here's some scilab code that does the cross-correlation and then finds the inter-sample peak by doing a simple (ish) parabolic interpolation.

The plot shows the absolute relative error between the truth and the estimate.

// 25684
vals= [];
taus = [-0.5:.005:0.5] + 10^-9; // Add a small amount to not have a zero.

for tau = taus

    N = 1000;
    t = 0:N-1;
    omega=0.1902834909;
    phi = 2*%pi*0.3823483924;

    x = sin(omega*t + phi);
    x_tau = sin(omega*(t-tau) + phi);

    xc = xcorr(x_tau,x);

    [mx,ix] = max(xc);

    x_vals = [-1 0 +1];
    y_vals = xc(ix + x_vals);

    //https://stackoverflow.com/a/717791/12570
    den = (x_vals(1) - x_vals(2))*(x_vals(1) - x_vals(3))*(x_vals(2) - x_vals(3));
    A = x_vals(3)*(y_vals(2) - y_vals(1)) + x_vals(2)*(y_vals(1) - y_vals(3)) + x_vals(1)*(y_vals(3) - y_vals(2))
    A = A/den;

    B = x_vals(3)^2*(y_vals(1) - y_vals(2)) + x_vals(2)^2*(y_vals(3) - y_vals(1)) + x_vals(1)^2*(y_vals(2) - y_vals(3));
    B = B/den;

    C = x_vals(2) * x_vals(3) * (x_vals(2) - x_vals(3))* y_vals(1) + x_vals(3) * x_vals(1) * (x_vals(3) - x_vals(1))* y_vals(2) + x_vals(1) * x_vals(2) * (x_vals(1) - x_vals(2))* y_vals(3);
    C = C/den;

    x_v = -B / (2*A);
    y_v = C - B*B / (4*A);

    vals= [vals; tau x_v y_v];

end

clf
plot(taus,abs((vals(:,2)' - taus)./taus));
xlabel('True value of delay');
title('Error in using parabola peak as subsample delay estimate');


Getting the bitstream of an unknown RF signal I am trying to get the bitstream of a signal sent over 433MHz. I have installed SDR# (I'm on Windows), connected my NESDR SMART, and I identified the signal here:

I can see there is only one peak with progressively decreasing peaks, so I assume this is AM. Then I set the bandwidth to the width of the central lobe. I finally made a recording of the "Audio" and opened the WAV into Audacity but I could not make sense of the sines (I was expecting a digital stream, since the modulation was specified) so here is the I/Q data instead.

This is an extract of the body of this waveform:

I have also installed DSD+ and the SDR# plugin as well as the virtual cable between the two, but after setting up SDR# the waveform does not budge so nothing comes up.

What sort of encoding is this? How can I read the bitstream? I am open to other things to try to get to the bottom of this.

First I found the sine wave extrema (black markers at the bottom) and stretched and superimposed a sine wave graph (red) on top of the signal graph with the extrema in the same locations. The signal matches (0) the sine wave over some of the non-zero segments, and over the other non-zero segments it matches the sine wave with a sign flip (1). This gives a binary sequence.

A sign flip is equivalent to a 180 degree phase shift, so this is a form of binary phase shift keying (BPSK, 2PSK, also phase reversal keying, PRK) with the useful symbol interval duration (of each sinusoid piece) seemingly equal to the guard interval duration (of each zero-valued segment that separates symbols). It is unusual that here the symbol rate is higher than the sine wave frequency, with a ratio of about $20/7\approx2.86$.


Improving an auto-correlation based guitar pitch detector I've seen many questions on this forum regarding pitch detection for musical instruments (commonly guitar), and spent a while reading through the answers to create a basic implementation of auto-correlation to make an Android guitar tuner.

This is the algorithm I'm using (implemented in Java on an Android phone):
1) Record a short[4096] array of raw audio data from an Android phone's microphone at 48kHz sample rate
2) Apply a Hanning window to the raw audio data
3) Zero-pad the result to double the length (8192)
4) Apply auto-correlation with FFTs
5) Normalize the auto-correlation result
6) Get the periodicity from the peak bin indexes

My problem is that with an actual guitar it is not robust (around 50% accurate at best), and I don't know how to filter noise either (without any loud noise, just ambient white noise, it outputs garbage frequencies).

It fares much better when I whistle at it, or play a generated sine tone from a computer, but that's expected.

In my search for ways to improve my implementation, many answers (on this forum and others) usually point to using a better algorithm like YIN, but I think its more likely that I've made a mistake in my implementation.

Does anything seem obviously incorrect from the algorithm I posted? Are there tweaks I can apply?

1) Increase the Android microphone sensitivty
2) Band-pass the raw data before the algorithm
3) Record more samples
4) ??? And does anybody have any idea how I could filter out noise by using the auto-correlation result?

Thanks in advance.

I realized I had forgotten to make an update to this.

I ended up following @robert bristow-johnson's recommendations in the comments.

I used time-domain McLeod Pitch Method in a live application (records audio in a loop) which works very well.

You can see the source code here:
https://github.com/sevagh/Pitcha


Find large distinctions in Video Assume you have 77777 pictures that are mostly noise.
They are about magnetic anisotropy.

There is one picture among these pictures which has high contrast difference between two large areas (some sort of magnetic gap because of magnetic anisotropy).
One area is gray and the other area is dark.

The pictures are mostly like this (grey or dark)

and

There are some dirty things in some pictures which are just rubbish on the lens.
So only big differences are essential.
Example what the contrast difference should look like

These pictures are about magnetic anisotrophy with Fe and Boron.

How can you find the distinction among these 77777 pictures?
I have went through all of these pictures but my method is errorprone, since I just browse them through.

How can you assist browsing the find the big and large contrast difference in pictures?

Any tool is ok for me: Matlab, Python, ..., pseudo-code too.

Implement a low pass 2D filter that will reduce greatly the noise.
Following an example that you could complete yourself.
You could subtract the filter bias to keep an average of 0 "addition" by the filter, controlling the dynamic range.

Example of filter

         1/16  1/8   1/4  1/8   1/16

    1/16  1/8   1/4   1/2  1/4   1/8     1/16

1/16  1/8   1/4   1/2   1    1/2   1/4     1/8     1/16

1/8   1/4   1/2   1     2    1     1/2     1/4     1/8

1/16  1/8   1/4   1/2   1    1/2   1/4     1/8     1/16

   1/16  1/8   1/4   1/2  1/4   1/8     1/16

        1/16  1/8   1/4  1/8   1/16

You might be able to clearly see the desired image.

Than run an Edge a 2D edge detector which will enhance the location of the step.
You could try running two 1D filters: -1 1 in horizontal and vertical direction. If you know the relative "black"/"white" areas you could select which to use of the optional -1 1 or 1 -1.

You might be able to use frequency domain filtering for the low pass 2D filter.
If you are not familiar with this type of challenges and this is not exercise you need to get a consultant to help you or study the subject in more depth.

Depending of the contrast of the "step" you might be able to detect automatically the "line" created by the edge detector.


Does a closed loop system with a reference of its previous iteration always reach a point in which it adventually loops, even if that loop is stable? For example you have a $F(\cdot)$ which starts off taking some input $a$, which describes some set of objects and their states, and creating a output $F(a)$, which can be renamed $F_1$, that also describes a set of objects and their states, which is then input into $F$, to get $F_2 = F(F_1)$, or $F(F(a))$.

If this is repeated for an extensive period is it guaranteed that that the values of the individual iterations will be loop? For example $F(F_{x+n})=F(F_{x+2n})$ etc and $F(F_{x+m})=F(F_{x+2m})$ where $x$ is the start of the cycle and $m$ and $n$ are the frequency of the cycle.

As is usual with these sort of questions, the answer is "it depends".

For example, suppose
$$
F(x) = x
$$
then the loop is immediate.

For example,  suppose
$$
F(x) = x/2
$$
then the loop never occurs unless $x = 0$.

The same will be true for more complex $F$s.

For example
$$
F(x) = \mathbf{A} x
$$
where $\mathbf{A}$ is a matrix. If $\mathbf{A}$ is a rotation matrix of a certain type (so that it rotates by $360/n$ degrees where $n$ is an integer) then the period will be $n$. If $\mathbf{A}$ is something else, then all bets are off.


Simulink for Amplitude Modulation I am using Matlab Simulink to analyse Amplitude Modulation in Noise. I attached my system below.

My problems:

1)In the mixer I am using a bandpass filter and another one after the envelope detector and I need to determine the minimum and maximum frequencies allowed depending on these parameters( I attached screenshots too) since I tried using mathwork site but it was not helpful regarding these values.

2)Also I need to know what needs to be changed in the model if the AM channel has a different carrier frequency.

3)What does 999e3 evaluate to?

Thanks a lot!

1) Do I fail to understand the question? The allowed frequency range of a bandpass filter goes from Fpass1 to Fpass2. The range between Fstop1 and Fpass1 and from Fpass2 to Fstop2 is the transition range where the filter goes from fully attenuated to fully "open" so there shouldn't be any strong interferers in this range.

2) Your band select filter before the down-mixer. It usually has your carrier frequency in the middle of Fpass1 and Fpass2 and is as wide as your expected signal is.

3) 999000


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

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

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

My guess is that somewhere in the processing you are using an AR (or ARMA) model that is being over / under fitted and that is what is causing the issue.

If I make a fake signal like this in R:

group_delay <- 84
T <- 350

omega <- 2*pi*0.252352890
tau <- 0.005
phi <- 2*pi*0.0989038
t <- 1:T

bruker_signal <- rep(0,T)
t_index <- seq(group_delay, T)
bruker_signal[t_index] <- exp(-(t_index-group_delay)*tau)*(
    sin(omega*t_index) + sin(1.5*omega*t_index + phi)) +
     rnorm(length(t_index),0, 0.1)

i.e. it's a noisy, delayed, damped sum of two sinusoids and then do AR modelling of the spectrum with different orders, then I get something like the smiles / frowns you are seeing.

An earlier attempt I made used a low-pass filtered version of the signal added to the same thing, modulated to $f_s/2$.

smile_filter <- butter(3,0.1)
smile_signal <- filter(smile_filter, bruker_signal)
plot(abs(fft(bruker_signal+2*smile_signal + 2*smile_signal*rep(c(1,-1),     lenght.out=T)))[1:175], type='l')
title('Attempt at Smile Spectrum ')

but that seemed way noisier than what you're after.


Extract breaths from respiratory signal, i.e. peak detection I am trying to extract breaths from a respiratory signal, which is essentially peak detection. See the example below for what such a signal looks like. This particular signal looks rather clean, however, lower amplitude noise is very common (both high and low frequencies).

I have tried the ones mentioned here, however, those seem to require that the height of the peak can be described by its y-value. As you can see below, the signal is very prone to baseline wander, and the top of a peak can very well fall below 0.

Any suggestions for a peak detection algorithm?

I have a slick trick which may do what you want.

Do an exponential smoothing in the forward direction, call it F.

Do an exponential smoothing in the backard direction, call it B.

Take the difference of these two D = (B - F) / 2.  (Rescaling by half isn't necessary)

This will do two things for you:

1) Smooth out the noise

2) Convert peaks (and troughs) to zero crossings, no matter what their height.

You can find details in my blog article "Exponential Smoothing with a Wrinkle".  There is a "ramp up" and "ramp down" zone at the beginning and the end that you should not use.  They are fairly short and their length depends on your smoothing factor.

Hope this helps,

Ced


Clock recovery using Mueller and Muller adds noise affecting EVM or SNR (Two cases - GNU Radio & Python code) The implementation of the Mueller and Muller clock recovery that I am seeing adds a large amount of 'noise' or degradation to the Error Vector Magnitude (EVM) to the symbols on the constellation. This is done independently in GNU Radio and a manual script in Python - giving similar results. Is this normal?
I understand that I need to adjust the parameters in this algorithm but I cant get anywhere near an optimisation that is adequate, its too noisy.
In both implementations GNU radio and python , the only impairment to the original data signal is the performance of the RRC pulse shaping filters. You can see that's its an almost perfect received constellation with an EVM < -60 dB using a filter taps of 300 and roll off factor of 0.35. See below image.

However, once I add the M&M algorithm I have no luck getting anywhere close to this and it generates a lot of noise. See below image.

It somewhat makes it difficult to assess the true extent of impairments I may later add in the model, like AWGN. Because if I add AWGN to this then the EVM or SNR just gets worse when the M&M is being used.
Here is the GNU radio flow.
EDIT Note that in my python code the Rx RRC filter is performed before the M&M, but in GNU radio I do it before an RRC because having an RRC before the MM block doesn't work give any result in the flow below.

Below is the M&M python algorithm, its from pySDR.com
samples = input_signal
interp = 16
samples_interpolated = signal.resample_poly(samples,interp,1)
mu = 2 #2 # initial estimate of phase of sample
out = np.zeros(len(samples) + 10, dtype=np.complex)
out_rail = np.zeros(len(samples) + 10, dtype=np.complex) # stores values, each iteration we need the previous 2 values plus current value
i_in = 0 # input samples index
i_out = 2 # output index (let first two outputs be 0)
while i_out < len(samples) and i_in < len(samples):
    out[i_out] = samples_interpolated[i_in*interp + int(mu*interp)] # grab what we think is the "best" sample
    out_rail[i_out] = int(np.real(out[i_out]) > 0) + 1j*int(np.imag(out[i_out]) > 0)
    z = (out_rail[i_out] - out_rail[i_out-2]) * np.conj(out[i_out-1])
    zz = (out[i_out] - out[i_out-2]) * np.conj(out_rail[i_out-1])
    mm_val = np.real(zz - z)
    mu += sps + 0.3*mm_val
    i_in += int(np.floor(mu)) # round down to nearest int since we are using it as an index
    mu = mu - np.floor(mu) # remove the integer part of mu
    i_out += 1 # increment output index
out = out[2:i_out-1] # remove the first two, and anything after i_out (that was never filled out)

Unlike the Gardner Loop, the M&M synchronizer should be performed after the RRC filter in the receiver for best performance.  With cases of high RRC alpha, the M&M won't work as expected without the complete Raised-Cosine filtering (RRC in transmitter followed by RRC in receiver) as the slope of the error term will reverse, with high self-noise, as I demonstrate in this plot below:

The blue line is the discriminator output showing average error versus offset, and the red line is the standard deviation of that error showing the self-noise. In the first case the slope of error output versus offset actually inverts in proximity of the lock point (so it would work, barely, if the feedback is inverted but with very high self-noise and low slope so very low timing SNR otherwise it would lock interestingly at 1/4 sample offset).
Specific to the OP's case, below is the same chart with a lower (and typical) RRC alpha where the discriminator curve is quite similar but the self-noise is much higher at the lock condition if the M&M is used prior to the second RRC filter.

Both the Gardner Loop and M&M will have self noise, but the self noise is typically averaged out in the recovery loop to be well below the receiver noise floor once the accumulated time sample location is determined for each new sample.
Why the difference between Gardner and M&M?
It is interesting that the Gardner performs better using the samples prior to the second RRC filter while the M&M synchronizer is better after. This is because the M&M synchronizer operates on one sample per symbol using sample locations at the decision locations and thus depends on the slope of the waveform in vicinity of these locations to determine timing error. The Gardner operates using 2 samples per symbol and operates by approximating the slope at the zero crossings, which have less timing jitter before the second RRC filter as we see in the eye diagrams below (notice the wider range of zero-crossing transition times in the waveform after the second RRC filter):

Generally for timing detection our best estimates should be obtained from the zero crossings where the slope and therefore sensitivity to timing error would be highest. This is true on average and is one of the reasons the Gardner Loop has better SNR performance in low SNR conditions (and M&M has better performance in high SNR conditions since it uses decisions in its computation, which eliminate noise if the decisions are correct). However, the competing requirement of having no inter-symbol interference at the data decision location results in increased zero-crossing jitter at the waveform locations of interest for timing detection. One solution to this if self-noise was an issue (it typically isn't as the timing loop bandwidth is so much smaller as to average it out) is the use of pre-filters which can equalize the ISI at the zero crossing locations for timing detection, on a separate data-path from the data which has zero ISI at the symbol detection locations as desired for best data-recovery performance.


Analogous filter to Kalman that maximizes mode (as opposed to minimizing variance) I may have a potential application where maximizing the mode (as opposed to typically minimizing the variance) would be useful for state estimates. The situation may arise from skewed distributions where the mode is of more interest than the mean. I was wondering if anyone has seen a filter like this before? Thank you.

Josh

You can use a recursive Bayes filtering (RBF) approach to do maximum a posteriori (MAP) estimation. While MMSE estimators look at the conditional mean, MAP estimators look at the conditional mode. See these notes for some details. In order to do RBF, you should look into particle filters / sequential monte carlo methods. There are plenty of books on this topic as well as review papers.

Alternatively, if you assume that your Kalman filter has white Gaussian noise independent for the process and measurement, and an independent Gaussian initial state, you can prove that the MAP estimator coincides with the Kalman filter. [ A Kalman filter is a RBF for the nice Gaussian case. ]

A good paper to look at is S. Godsill, A. Doucet, M. West's “Maximum a posteriori sequence estimation using Monte Carlo particle ﬁlters,” Ann. Inst. Stat. Math., vol. 52, no. 1, pp. 82–96, 2001.

The actual article is behind a paywall, but certainly, many articles on related topics have been written (You can just google Maximum A Posteriori estimation with particle filtering or recursive bayes filtering for details).


OpenCV - What type of a result does cv2.goodFeaturesToTrack() return? Dear citizens of Stackexchange,

I am a 3rd grade Electrical and Electronics Engineer student, who is working on an image processing project using OpenCV 3.2 and Python 2.7.

What we would like to accomplish is to build up a system which can compensate the camera using Lucas Kanade's Optical flow method.

On the sample that was presented on the OpenCV's website, one of the inputs of the function calcOpticalFlowPyrLK() was cv2.goodFeaturesToTrack.

I wanted to change that input to an output I have acquired from my another function. But the outputs of those two functions are mainly different in type.

What I would like to ask is, what type of topology is the following output format (matrix, array, etc.) that we get from cv2.goodFeaturesToTrack()?:

[[[ 582.  408.]]
 [[ 453.  373.]]
 [[ 483.  395.]]
 [[ 447.  298.]]
 [[ 589.  371.]]
 [[ 436.  358.]]
 [[ 624.  370.]]
 [[ 624.  413.]]
 [[ 470.  395.]]
 [[ 635.  378.]]
 [[ 630.  417.]]
 [[ 498.  395.]]
 [[ 427.  354.]]
 [[   6.  325.]]
 [[ 566.  381.]]
 [[  20.  323.]]
 [[ 532.  400.]]
 [[ 568.  374.]]
 [[ 457.  321.]]
 [[ 460.  413.]]
 [[ 471.  317.]]
 [[ 524.  399.]]
 [[ 586.  417.]]
 [[ 251.  184.]]
 [[ 289.  428.]]
 [[ 435.  342.]]
 [[ 635.  393.]]
 [[ 590.  387.]]
 [[ 160.   20.]]
 [[ 431.  296.]]
 [[ 449.  324.]]
 [[ 439.  371.]]
 [[ 466.  372.]]
 [[ 481.  465.]]]

The cv2.goodFeaturesToTrack() is mainly used for detection of corners. See THIS PAGE

If you observe the type of the value returned by  cv2.goodFeaturesToTrack(), it is numpy.ndarray, which is a multi-dimensional array according to THIS DOCUMENT

So the function cv2.goodFeaturesToTrack() returns all the locations of the corners present in the gray scale image.

This is what a typical output would look like, the corners are marked in blue:


FIR filter with linear phase, 4 types I know there are 4 types of FIR filters with linear phase, i.e. constant group delay: (M = length of impulse response)

Impulse response symmetrical, M = odd
Imp. resp. symmetrical, M = even
Imp. resp. anti-symmetrical, M = odd
Imp. resp. anti-symmetrical, M = even

each with its traits. Which of these types is most commonly used in FIR filter with linear phase design and why? :)

When choosing one of these 4 types of linear phase filters there are mainly 3 things to consider:

constraints on the zeros of $H(z)$ at $z=1$ and $z=-1$
integer/non-integer group delay
phase shift (apart from the linear phase)

For type I filters (odd number of taps, even symmetry) there are no constraints on the zeros at $z=1$ and $z=-1$, the phase shift is zero (apart from the linear phase), and the group delay is an integer value.

Type II filters (even number of taps, even symmetry) always have a zero at $z=-1$ (i.e., half the sampling frequency), they have a zero phase shift, and they have a non-integer group delay.

Type III filters (odd number of taps, odd symmetry) always have zeros at $z=1$ and $z=-1$ (i.e. at $f=0$ and $f=f_s/2$), they have a 90 degrees phase shift, and an integer group delay.

Type IV filters (even number of taps, odd symmetry) always have a zero at $z=1$, a phase shift of 90 degrees, and a non-integer group delay.

This implies (among other things) the following:

Type I filters are pretty universal, but they cannot be used whenever a 90 degrees phase shift is necessary, e.g. for differentiators or Hilbert transformers.
Type II filters would normally not be used for high pass or band stop filters, due to the zero at $z=-1$, i.e. at $f=f_s/2$. Neither can they be used for applications where a 90 degrees phase shift is necessary.
Type III filters cannot be used for standard frequency selective filters because in these cases the 90 degrees phase shift is usually undesirable. For Hilbert transformers, type III filters have a relatively bad magnitude approximation at very low and very high frequencies due to the zeros at $z=1$ and $z=-1$. On the other hand, a type III Hilbert transformer can be implemented more efficiently than a type IV Hilbert transformer because in this case every other tap is zero.
Type IV filters cannot be used for standard frequency selective filters, for the same reasons as type III filters. They are well suited for differentiators and Hilbert transformers, and their magnitude approximation is usually better because, unlike type III filters, they have no zero at $z=-1$.
In some applications an integer group delay is desirable. In these cases type I or type III filters are preferred.


Software implemented PLLs for signal synchronization, why use them? In short I'm recreating a program that reads 3 sinusoidal signals through an ADC performs some manipulations and then reproduces these signals through a DAC. In order to reproduce these signals I need the instantaneous phase of one of these manipulations. So I recreated a charge pump zero crossing pll in software.

I don't understand why you would go to the trouble of working out the analogue loop filter transfer function, convert it to the z domain, when in software when I could just implement a program that reads the amount of samples between zeros, and the instantaneous phase is just.

(current sample#-last zero sample#)/(pi*(number of samples between zeroes))

Now I've obviously missed something, but I don't know what to google to find out what.

For very high SNR signals locating the zero crossings will work pretty well.  You use PLL's when the signal doesn't have a high SNR, like the following.

As you can see by inspection, a zero-crossing algorithm wouldn't have a prayer of working on this signal.  A PLL, on the other hand, could do just fine.  That is, by the way, a tone with an SNR of 8 dB.


Linear prediction and filter stability I am currently trying to implement an iterative block-wise algorithm in which AR coefficients are computed for each block.
There is a issue with my code where values get too large, and I was wondering how linear prediction handles stability issues. For an onset block for example, the coefficients for a low order IIR-filter may result in an unstable filter, at least from my intuition.

My question is:

Is the linear prediction method able to produce instable sets of autoregressive filter coefficients?

I realize this might be a stupid question, and I feel like I missed something.
Please note that I would like to hear about linear prediction, I am not asking for someone to fix my code problem (yet).

Ok, I found out from a book that linear prediction always results in minimum phase filters, which implies stability.

This can be shown using the mean-squared error (MSE) criterion which is to be minimized in the autocorrelation method. The proof is done by contradiction; if an all-pole filter had a single pole outside of the unit circle, the derivative of the MSE in regards to the pole radius is always greater than zero. In terms of the MSE, this implies the pole position is not optimal, so it cannot be a result from linear prediction (under ideal conditions, I assume).

The actual proof can be found here.


Z transform of a function with delay? i have this open loop system , and i've been asked to find out the response $C(kT)$ due to a unit step input.

I am able to find the transfer function without the delay unit i.e $$\frac{C(z)}{R(z)}=\left(1-z^{-1}\right)\cdot \mathcal{Z}\left\{\frac{1}{s(s+a)}\right\}$$
i can find this ZT easily and find $C(kT)$

But how to deal with the exponential term i.e $e^{-s}$ ? then my pulse transfer function becomes
$$\frac{C(z)}{R(z)}=\left(1-z^{-1}\right)\cdot \mathcal{Z}\left\{\frac{e^{-s}}{s(s+a)}\right\}$$
Can i also consider it as a delay ? and take it out of the braces ? i.e
$$\frac{C(z)}{R(z)}=\left(1-z^{-1}\right)\cdot z^{-1} \cdot \mathcal{Z}\left\{\frac{1}{s(s+a)}\right\}$$

In my opinion, the most straightforward way to solve this problem is to consider the time domain signal as it passes through the signal chain. A unit step at the input when sampled and put through a ZOH becomes just another continuous-time unit step. Consequently, the Laplace transform of the continuous-time signal just before the sampler at the output is given by

$$C(s)=\frac{e^{-s}}{s(s+a)}=\frac{e^{-s}}{a}\left[\frac{1}{s}-\frac{1}{s+a}\right]\tag{1}$$

which corresponds to the time domain signal

$$c(t)=\frac{1}{a}\left[1-e^{-a(t-1)}\right]u(t-1)\tag{2}$$

Sampling $c(t)$ at $t=kT$ gives

$$c(kT)=\frac{1}{a}\left[1-e^{-a(kT-1)}\right]u(kT-1)\tag{3}$$


Object Classification using pyramid match kernel I have been reading up on the pyramid match kernel as an alternative to using the bag-of-words model for object classification using a SVM. The bag-of-words model provides a model for transforming a set of features to a single feature vector, however im confused as to how the pyramid match kernel can be used for this task, for instant if a SVM has been trained using a pyramid match kernel how do I present a test image represented by a set of feature vectors to the SVM as, as far as I can tell the kernel generates some feature space F(x) but provides no mapping to this space. Could someone please give a brief description of how this works?

The pyramid match kernel does not operate on the sets of feature vectors directly.  It operates on multi-level histograms of feature vectors. Let's say you are using the SIFT descriptors, which live in 128-dimensional space.  First you divide each dimension into 2 bins, which divides the entire feature space into $2^{128}$ bins. Then you count how many features fall into each bin.  This is the first level of your histogram.  Then you divide each dimension into 4 bins, and so on. Luckily, this histogram is going to be very sparse, so you only need to allocate the memory for non-zero bins.

You do this for the training images, and the images that you want to classify.  Then the pyramid match kernel is simply the intersection of two histograms.


Using MIMO with DS-CDMA and spread spectrum I have a question about using MIMO with DS-CDMA. suppose we are using 4 TX antennas, and 4 RX antennas, then the Reyleigh channel generated randomly  of MIMO system will be of dimension  H = [4x4],  assuming also the transmitted signal using QAM modulation is x done as below:

x =
[0.7 + 0.7i;
0.7 - 0.7i;
-0.7 + 0.7i;
-0.7 -0.7i];

So the received signal r is supposed to be r = H*x , which is of dimension of [4x1].

Now, suppose we are spreading that signal with PN/walsh code of dimension [4x1], by using the function of "kron" in MATLAB, so the dimension of the new transmitted signal will be [16x1], equivalent that x_1 = kron(x,c); where c is the code used to spread the signal before transmission.

My question is, how can we transmit that signal after spreading it x_1 over the noted channel H of dimension H with dimension [4x4]?

Thank you

Assuming one antenna transmits one symbol per time unit, then 16 symbols require 4 time units to be out. Then it is simply that

r_1 = H_1 * x_1(1:4)
r_2 = H_2 * x_1(5:8)
r_3 = H_3 * x_1(9:12)
r_4 = H_4 * x_1(13:16)

If channel H is fixed during these 4 time units,

[r_1 r_2 r_3 r_4] = H * [x_1(1:4) x_1(5:8) x_1(9:12) x_1(13:16)];

or

r = reshape(H*reshape(x_1,4,[]),[],1);

This can be generalized to any MIMO size and to any spreading size.


How to express STFT and ISTFT as a 1d convolution and 1d deconvolution in tensorflow/keras I'm trying to implement this paper in tensorflow and keras. At the end of section 3 it says.

Note that although weighted SDR loss is a time-domain loss function,
it can be backpropagated through our framework.
Specifically, STFT and ISTFT operations are implemented as 1-D convolution
and deconvolution layers consisting of fixed filters initialized
with the discrete Fourier transform matrix.

How would I go about implementing an STFT and an ISTFT as a 1-D convolution and 1-D deconvolution respectively? For example how would the window size and the hop length correspond to parameters of a convolution?

Conversely, could I just use tensorflow's built in STFT and inverse_STFT ops instead?

I found an implementation of STFT based on conv1d in pytorch here:

https://github.com/huyanxin/phasen/blob/master/model/conv_stft.py

edit: Actually, the phasen repository took the STFT code from https://github.com/pseeth/torch-stft
edit2: Asteroid has an alternative implementation of STFT and iSTFT: https://github.com/mpariente/asteroid/


What is special about the frequency $\omega_0=\pi$ that suddenly causes rate of oscillation decrease? In Alan Oppenheim's book Signals and Systems a comparison is made between the properties of discrete-time and continuous-time complex exponential signals in section 1.3 pg. 26.  Specifically it says:

  The continuous-time complex exponential $e^{j\omega_0t}$ has two
  properties:
      1.  the larger the magnitude of $\omega_0$, the higher is the rate of oscillation in the signal
      2. $e^{j\omega_0t}$ is periodic for any value of $\omega_0$

  Consider the discrete-time signal $e^{j(\omega_o+2\pi)n}= e^{j2\pi n}
e^{j\omega_0 n}=e^{j\omega_0 n}$. Because of the periodicity implied by this equation, the signal $e^{j\omega_0 n}$ does not have a continually increasing rate of oscillation as $\omega_0$ is increased in magnitude.  Rather, as illustrated in Figure 1.27, as we increase $\omega_0$ from 0, we obtain signals that oscillate more and more rapidly until we reach $\omega_0=\pi$.  As we continue to increase $\omega_0$, we decrease the rate of oscillation until we reach $\omega_0=2\pi$, which produces the same constant sequence as $\omega_0=0$

My confusion is with the part that is in bold.  For the continuous-time case, we know that

$\omega_0 = \frac{2\pi}{T_0}$ and it makes sense to me that decreasing $T_0$ (more frequent oscillations) would cause $\omega_0$ to increase.  Similarly, for the discrete-time case we know that

$\omega_0 N = 2\pi m$ for $N,m \in \mathbb{Z}$ Solving this equation for $\omega_0$, we get

$\omega_0=\frac{2\pi m}{N}$.  Just like in the continuous-time case, decreasing N should cause $\omega_0$ to increase.  So I am not seeing what is so special about $\omega_0=\pi$, the rate of oscillations suddenly start decreasing. I don't see that in the math if you decrease N.

As a consequence of the sampling theorem, the frequency spectrum of a discrete random variable is periodic, but it turns out that signals that are inherently discrete in time and not derived from an A/D converted continuous time signal are also periodic in frequency.  As a consequence, as the frequency is increased, it goes from the positive portion of one periodic interval, to the next, i.e. it transitions to the negative image of the next interval.

In the old days, DSP was mostly a sub discipline of Electrical Engineering, and DSP was focused on doing things digitally that were previously done analog.  In order to talk about DSP, without having to reference, a sampling frequency, a convention called unit sampling was commonly used, and a period of the spectrum is from -1/2 to 1/2.

As time progressed it was noted that DSP wasn't directly tied to processing analog derived signals and unit sampling was somewhat artificial, so it was recognized that $2\pi$ was a natural periodic interval, an $2\pi$ fits between $-\pi$ and $\pi$, which is why $\omega=\pi $ is a transition point in the spectrum.


Separating signal from noise in similar frequency range As shown in the figures, both of my signal and noise have the same frequency domain, how can I separate them?

As you have noticed, it is quite unlikely that you can filter your data with a standard frequency filter, thus a linear filter.

Unless the signal and the noise are well separated in time, unless you have a model, you will probably need to resort to other techniques, such as nonlinear filtering, time-frequency or time-scale filtering, or add additional information or assumptions, such as thought signal or noise properties (to name a few: gaussianity, positivity, independence, sparsity), and bend to blind source separation, latent variable analysis, nonlinear matrix separation, etc.


High pass filter not attenuating signal at start of input I have implemented digital high pass filter with below specs:
Sampling frequency = 5333 Hz
Cutoff frequency = 10 Hz
Order = 3

We are facing problem like below:
As per result shown above signal gets attenuation after $1000$ data points. How can solve this kind of issue?
Any help/ suggestion would be great and appreciable! Thanks in advance.

Any filter with such a low cutoff (as compared to the sample rate) will lots of time domaing ringing.
The best way of adressing this is to properly initialize the state of the filter using known properties of your system & signal, but that's fairly tricky.
In your case, there may be a shit cut though: I'm guessing that you have a large bias or DC offset in your original signal. Just subtract the mean (or maybe the mean over the first samples or so, if the bias is drifting). This should subtantially reduce the ringing.


Square law demodulator. How can we eliminate the term containing $m^2(t)$? We are trying to extract the message signal from the AM signal using square law demodulator.

In the final expression, $V_2(t)$, the term encircled in blue has to be retained, and the others have to be eliminated. The terms encircled in green can be eliminated by low pass filter(the career frequency being high). The constant term can be eliminated with the help of coupling capacitor as mentioned below the derivation. But the problem is that the spectrum of the term $\frac{k_2 {A_{c}}^{2}{k_{a}}^{2}m^2\left ( t \right )}{2}$ will lie around the origin. It is just that the width of the spectrum of this signal component will be twice of that of $m(t)$[Multiplication in time domain is convolution in frequency domain]. The spectrum of the the component $k_2{A_{c}}^{2}k_am\left ( t \right )$ will also lie around the origin. Will the spectrum of $k_2{A_{c}}^{2}k_am\left ( t \right )$ and  $\frac{k_2 {A_{c}}^{2}{k_{a}}^{2}m^2\left ( t \right )}{2}$ not interfare? How can then we extract the message signal?

$|m(t)| \leq 1.0$, so for analog signals, ${m^2(t)}$ is going to be small most of the time.  Also $k_a < 1.0$, so ${k_a}^2 \ll 1.0$.  You can ignore the term usually.  If you can't, then use something other than a square law demodulator.


CIC Filter vs Simple Downsampler I' trying to implement a QPSK modulation/demodulation.

At one point I have to decimate my signal by 4 for an easier implementation. I've read everywhere that the best way to do it is to use a CIC filter followed by a decimator. So I've read some more and I decided to use the following:

CIC Filter Length: 2
CIC Filter Order: 5
CIC Filter Decimator: 2
CIC Filter Stages: 2

So basically what I do is:

Filter the signal with a FIR filter and coefficients [1 5 10 10 5 1].
Downsample by 2.
Filter the signal with a FIR filter and coefficients [1 5 10 10 5 1].
Downsample by 2.

My question is:

  What is the benefits of a CIC Decimator compared to directly downsampling by 4 and thus disregarding 3 samples out of 4?

I simulated this in order to find what the difference in performance was. As you can see in the two figures below the CIC filter seems to have a strong degrading impact on the BER which makes me wonder why it is used.

If you are only decimating by four, then a CIC filter is not the way to go.  CIC filters have two advantages and one big disadvantage.  The advantages are that they don't require any multipliers, just (large) adders, and they can efficiently decimate by really large factors.  The no-multipliers advantage was a big deal years ago when multipliers were expensive and not always available, but they are cheap now.  And since you are only decimating by a factor of 4, the ability to decimate by large factors does not help you either.

The one big disadvantage of CIC filters is that their frequency response is not flat in the pass-band, so designers usually put a FIR filter after the CIC filter whose primary purpose is to compensate for the CIC filter's non-flatness.

My advice would be to continue decimating with one or two FIR filters, though I would recommend using a better filter than [1 5 10 10 5 1].  Try using Matlab's fdatool or Matlab/Octave's "fir1" command to design a better filter.


Does taking the Hilbert transform extract the envelope of audio signals? I want to extract the envelope of music signals with a sample rate of 44.1kHz.

I used the MATLAB command:  abs(hilbert(mysignals)), but the resulting signal is similar to the original signal, not the envelope.

Does taking the Hilbert transform extract the envelope of an audio signal? Or is my calculation incorrect?

That really depends on your definition of "enevelope" and what you need it for.

The Hilbert transform calculates the "analytic" signal, i.e. it calculates a matching imaginary part to a real signal by shifting the phase by 90 degrees in the frequency domain. It's reputation of calculating the "envelope" comes mainly from communication technology. It works very well with a narrow band signals, like amplitude or frequency or phase modulated signals. It's based on the simple fact that the analytic signal of a sine wave is a complex exponentional, the magnitude of which is 1. So the hilbert transform (or to be precise calculating the magnitude of the analytic signal) will eliminate the carrier quite well and leaves you with the modulation signal, or "envelope".

For broad band music, this is mostly not the case and the analytic signal isn't particularly helpful. A cymbal crash is basically a noise burst. The hilbert transform of this is just another noise burst and so is the magnitude of the analytic signal. You can't eliminate a carrier, if there isn't one to start with.

IF you are looking for something that roughly gives you "perceived loudness as a function of time", there many different ways to do this from simple lossy peak detectors all the way to sophisticated perceptual algorithms including cochlea models, binaural & temporal masking, nerve excitation models etc. It really depends on your application.


How to determine if a pixel is inside a contour? I have following picture

which I extract the contour from. The contour looks like following

I am wondering, for any pixel in the image, how do I determine if it is inside the contour or not? Say "inside" here means any pixel in the body.

If possible, I'd prefer simple, direct, and elegant algorithm. Thanks.

I am not in the image processing field, but here is an approach from my control loop background that I believe would apply as an elegant solution for determining if a pixel is inside a contour.  This would be to treat the image as a complex plane and then use the Cauchy Argument Principle (from 1831!!!) to determine if the pixel is inside or outside the contour.

The details and background are below, but as an algorithm you would do the following:

Treat your pixel of interest and your contour as a complex number, meaning a pixel at an x, y location is represented as x+ jy (I am an engineer so trained to you j instead of i).  Using complex number subtraction, subtract the pixel from every point on your contour (which is a new vector with magnitude and phase) and measure the phase of this new vector. Continue along the contour while keeping track of the minimum and maximum phase value in the results after the subtraction. Once the contour is complete (you have completely enclosed your pixel) determine the maximum - minimum phase. If it is 2$pi$ the pixel is inside your contour, if not then it is outside.

See figure below for an introduction to the Cauchy Argument prinicple, and then I will try to explain how this becomes the algorithmic solution above. From this you may define other similar algorithms depending on the details of your contour, such as you will notice that the phase in the example given where the contour phase keeps increasing will also keep increasing when the pixel is inside the coutour, but will increase and then decrease when it is outside.

To explain the principle I will need to explain (in simplest terms) the idea of poles and zeros and encirclement on the complex plane.  For your case we could simply identify the one pixel as a "zero" and proceed from there, which makes the explanation even simpler, so I will avoid any details on encircling poles.
(My plot is titled "Nyquist Plot Background" as this is the theory behind the Nyquist Plot that is used in classical control theory).

Given a ratio of polynomials, the poles are the values from a complex plane that would result in the ratio being infinity, and the zeros are the values that would result in the ratio being zero. Typically the polynomials are factors into first order functions to readily reveal the pole and zero locations. For example:

$$H(s) = \frac{(s+2)}{(s+3)(s+5)}$$ would have a zero at s= -2 and two poles at s=-3 and s=-5. These are the values that would make the ratio go to zero and infinity. At this point for a test of a single pixel, we can proceed with the idea of a single zero so will use that going forward.

In terms of an image with two axis x and y, x would correspond to the real numbers in this case, and y would correspond to imaginary numbers (the solution in the end will be very simple and elegant, but we need to use complex numbers to use it). So in the case above, the zero would be on the image at x= -2 and y = 0.  To position a zero at an arbitrary location with x and y values we would use $(s - (x+jy))$ as the zero.

That said, observe how the Cauchy Argument Principle and complex mapping works in the case of a single zero:

In the following figure, we have a mapping of a zero at location x= -2, y=0 defined by the function $$H(s) = s+2$$

The mapping of H(s) for the example contour shown is on the right. What the mapping is specifically is to take each point on the contour and subtract it from the zero. So for each point $(x_n, y_n)$ on the contour, and the pole at $(x_z, y_z)$ we would get the result $(x_n+j y_n) - (x_z+j y_z)$ which is a complex vector with the magnitude and phase just as shown by the blue vector on the left in the plot below. Note that the magnitude is simply the geometric length between these two points in 2D space, and the angle is the angle found by moving the zero to the origin; by subtracting the two, we are simply moving the zero to the origin which is the complete mapping in the case of the zero (simple!). The Cauchy Principle reduces to saying if you encircle the zero that is inside the contour, the mapping will encircle the origin in the same direction.

In contrast, if the zero is outside the contour, as in the plot below with the zero at x= -4, y =0, the mapping will not enclose the origin.

So using this, we can position your one pixel to be the zero as a complex number at the x,y position, and complete the mapping by subtracting each x,y position on your contour to each x,y position on your pixel. Now to algorithmically test that the result encircles the origin, note the phase for each of the two cases above. For the case of encirclement, the phase will go from 0 to $2\pi$ always. For the case when we don't encircle the origin the phase can never extend completely from 0 to $2\pi$.


Conceptually confused by LPC for speech: Do we synthesize by the inverse filter (FIR)? I've been following the process of using LPC to analyse the speech and then synthesize the speech by swapping out coefficients over time. I'm purely concerned with pitched vowels at this point.

lpc(speech_segment) gives us $a_{p}$ coefficients of an IIR filter, where $p$ is the order. This IIR filter has the frequency response of the speech vowel/formant present in speech_segment.

I have a pulse train to use as a source.

I see two ideologies for synthesis:

some texts talk about using the IIR filter to filter the source, which makes sense to me intuitively - juxtaposing the vowel's frequency response to my source's spectrum
Other texts talk about using the inverse of this IIR filter to form an FIR, and then filter the source somehow. I think I understand why this is invertible and stable (minimum-phase).

So approach 2 makes less sense to me from source-filter perspective, but it would be really useful because it means I can perform LPC analysis-synthesis with OLA convolution of audio blocks rather than implement the filter recursively.

Is my understanding of either of these two approaches, not quite right? If not, what exactly changes to make each possible.

Given speech samples, the LPC computation yields linear prediction coefficients $a_1,\ldots,a_p$. These describe a dependence model for a few samples. It is assumed that for short-enough collections of speech samples, the process is modeled well by an autoregressive stochastic process and that each sample is well-approximated by the following estimate based on the previous $p$ samples:
\begin{equation}
\widehat{x}[n] = a_1x[n-1] + \cdots + a_px[n-p].
\end{equation}
The error is then
\begin{equation}
e[n] = x[n] - a_1x[n-1] - \cdots - a_px[n-p].
\end{equation}
This "error filter" has coefficients $1, -a_1,\ldots,-a_p$.

For un-voiced speech, it turns out that discrete-time white noise is a decent model for the error signal. Keep this in mind.

Rather than send all the samples, communication systems send these linear prediction coefficients$^{\dagger}$. As mentioned in the other answer, to generate an approximation of the original speech, the synthesizer either generates a pulse train or a white noise realization and sends it though the filter that is the inverse of the error filter. This inverse is an IIR filter with frequency response
\begin{equation}
H(e^{i\omega}) = \frac{1}{1 - a_1e^{-i\omega} - \cdots - a_pe^{-ip\omega}}.
\end{equation}
The output of that IIR filter is an approximation of the original collection of speech samples.

$\dagger$ Vocoders do not really spit out linear prediction coefficients. They do not even spit out proxies such as reflection coefficients or line spectral pairs. They spit out integer indices of quantized versions of these coefficients. When these integers are received, the synthesizer grabs the appropriate entries in a lookup table and uses those as the coefficients.


Simple, streaming, lossless image compression Does anyone know of any image compression techniques with the following characteristics:

lossless
streaming - I want to compress on the fly, pixel-by-pixel.
low-memory overhead - I can afford to buffer a single line, but ideally not even that.
no dynamic dictionaries
"real-world" images only, so performance on "nasty cases" like chequerboards is not important
2-3x compression (5x-10x would be even better, but that's asking a lot I know)
can operate on 10-16 bit pixels (depending on my camera)

My images will be ~1k pixels wide, with pixel rates of ~20Mpix/sec.  The pixel depth will be something between 10 and 16 bits per pixel (depending on the choice of camera). Assume sub-16-bit pixel widths would be represented within a 16-bit word for now, rather than needing to be extracted from a continuous bit-stream.

Some form of delta+arithmetic coding perhaps?

You can consider using Huffyuv : http://neuron2.net/www.math.berkeley.edu/benrg/huffyuv.html

This is no great better than simple zip, but still slightly optimized for images.

Any Image related compression comes from the techniques like Vector quantization or Transform coding. In order to make use of transform such as DCT/Wavelet yet make it lossless you can think of JPEG-LS or JPEG2000 for compression. Only thing is, it is not streaming in your sense of definition.


What is the definition of small/large signal bandwidth? What is the definition of 'small signal bandwidth' and 'large signal bandwidth'?
If I google them, I only encounter them in an electrical context (transistors etc.). However, I am looking for a general definition (for any time signal).

In communications, whether a signal is narrowband (has a small bandwidth) or wideband (has a large bandwidth) is relative to the channel's impulse response.

If the channel behaves as a flat filter (constant gain) over the signal's bandwidth, then the signal is narrowband. If the channel introduces significant distortion (that is, if different signal frequencies see different gain) then it is said that the signal is wideband.


Digital lowpass filter causing aliasing I'm writing a biquad filter program in C++, using the LBJ Cookbook as the source for the formulas. https://webaudio.github.io/Audio-EQ-Cookbook/audio-eq-cookbook.html
The lowpass filter that I've made seems to half work. The frequencies above the cutoff are being attenuated, however, it is causing aliasing similar to what happens with frequencies above Nyquist's limit.
double output, a0, a1, a2, b0, b1, b2, w, alpha;
double pi = 3.14159265359;
w = (2 * pi) * (cutoff / sampleRate);
alpha = (sin(w)) / 2 * q;
a0 = 1 + alpha;
a1 = -2 * cos(w);
a2 = 1 - alpha;
b0 = (1 - cos(w)) / 2;
b1 = 1 - cos(w);
b2 = (1 - cos(w)) / 2;

output = (b0 / a0) * inSamples[0] + (b1 / a0) * inSamples[1] + (b2 / a0) * inSamples[2]
  - (a1 / a0) * outSamples[0] - (a2 / a0) * outSamples[1];

I'm not sure where my error is, so any help is greatly appreciated!

your output is definitely a linear combination of previous in- and output samples.
That means it's a linear time-invariant operation.
That means it can't introduce aliases.
What you see is something else than an alias, and I'll venture the guess you're interpreting sidelobes as aliases.


SineWave generation and BandPass filter : Audio Can any one suggest me a Bandpass filter in java which I would like to apply on generated sine wave of a particular frequency

As someone suggested in the comments, filtering a sine wave is just a change in amplitude and possibly phase, although you will have some edge effects at the beginning. I assume this is not simply an excersize in changing amplitudes and that you really need a real filter.

I don't know of an open source library in Java, although this book contains some simple source for a few filters (not sure it has a bandpass, but lowpass + highpass = bandpass).

I don't recommend the suggestion in the comments of using an FFT, although if you can perform a transform on your entire input that is an option, and one that many people find very intuitive. Keep in mind it's an O(n log n) operation, so it may not be feasible for large or even moderately-sized audio files. I explain here more on why EQ is usually not done in the frequency domain. Some people design their filters in the frequency domain and use the IFFT to create a filter, but this results in a design that is inferior in many ways. There are some rules of thumb that make this better, but there are better techniques, so this is not a well developed art.

You haven't specified your filter in any meaningful way, but assuming a basic biquad will suffice, checkout RBJ's audio EQ cookbook. If you need help implementing those filters, I've written a tutorial.

Finally, someone else on StackOverflow just asked about a bandpass filter in Java, and offered some not yet working code. Their design is based on a blackman windowed sinc filter. This is a high-quality FIR filter.


What does it mean by DN in salt-pepper noise? What is DN on page 245 of these notes shown in the slide below?

From the same course, you can find the following indication in an earlier set of slides:

  Input quantized image pixel values (integers): Digital Number (DN)

From this you can infer that

  Salt noise: DN = maximum possible

indicates that salt noise shows up as pixels with the largest possible pixel values (e.g. 255 if the picture is a grey scale picture represented by pixel in the [0-255] range), and that

  Pepper noise: DN = minimum possible

indicates that pepper noise shows up as pixels with the smallest possible pixel values (e.g. 0 if the picture is a grey scale picture represented by pixel in the [0-255] range).


What are the advantages of using fft2 over fft Assume we have a matrix x of size (8,8), . As known,  FFT(x) performs 1D-FFT transformation, column wise. However,  FFT2(x), performs 2D-FFT transformation.

In that case, what's the advantage of using 2D-FFT over 1D-FFT ?

They perform two different mathematical operations

FFT executes a 1-dimensional Discrete Fourier Transform one each column of the input matrix (or the first non-singleton dimension)
FFT2 executes a full 2-dimensional Discrete Fourier Transfrom on the entire matrix.

Which one you want to use depends on your specific application. One is not inherently "better" than the other, they are just two different things. In particular

$$\mathcal{F}_{2D}(X) = \mathcal{F}(\mathcal{F}(X)')'$$


What's the Difference Between UQI and SSIM Measures for Image Similarity / Quality? I'm studying image quality assessment (IQA) problem and I'm currently learning about quality measures UQI and SSIM, I know UQI is a predecessor of SSIM, but they look exactly the same for me...

The developers of both are the same hence the similarity is indeed "By Design".

The only difference is the addition of 2 constants in SSIM (C1 and C2).

The UQI:

The SSIM:

As the writers write in the SSIM paper:

Namely, UQI is a private case of SSIM for C1 = 0 and C2 = 0.


How to design a narrow band pass filter without DSP toolbox I've two signals X and Y. These are proximity sensors reading measuring vibrations of a rotating shaft. Now I want to plot orbit plots for these signals. I can plot unfiltered orbits with just X and Y. But since there is lot of noise in the signal I filtered it. Now there are dominant 1X frequency component and may be other sub frequency components. I want to plot the filtered orbit of 1X say 656 Hz. The sampling frequency I used is 1000000 Hz. I know band pass filter could give me only the 1X component. But I don't have DSP toolbox. Could someone help me how to design and implement narrow band pass filter without toolbox.
Any suggestion or hints or idea would be helpful.
Thanks in advance.

The steps that I would follow to handcraft a bandpass filter are these:

Figure out how wide (in Hz) you want the filter to be.
Create a low-pass filter with half the one-sided bandwidth as your desired filter bandwidth.  You can calculate the sinc function with basic functions, since $sinc(x) = \frac{sin(\pi x)}{\pi x}$.  You can adjust the bandwidth by adjusting the time scale in your sinc function.  If n is the set of all integers between -N and N, inclusive, and you make your filter taps equal to $sinc(bn)$, then the filter's normalized bandwidth will be equal to $b$.
Modulate the lowpass filter to the center frequency you want for your bandpass filter.  You do that by multiplying it, element by element, with a sine wave of that frequency.


MUSIC algorithm for frequency and DOA estimation $\lambda$ denotes the wavelength of the incident sinusoid, $\theta$ the direction of arrival (DOA), and x the distance between the sensor array elements. The signal is $x_{1}(n)....x_{k}(n)$
How do I proof/explain that estimating the DOA($\theta$) is the same as estimating the frequency of the signal (where $x_{1}(n)=x_{1}(m)$ instead where m is the time instant.

I understand that the MUSIC algorithm can be used to estimate the frequency of a signal, and the MUSIC algorithm can be used to estimate the DOA of a signal as well. They are both related as they are basically finding the peak of the power spectrum, but I'm not entirely sure how to explain this.

The estimation of frequency for one dimensional signal (dimension is time) is same as one dimensional space data for narrow band signal. To make it clear, consider you have n time samples (Nyquist criteria is met) data. The problem of frequency estimation from this n sample is same as the following case. Consider you have a n element linear array (equal spacing) and you have received a time stamp data from a narrow band signal  i.e you have n samples collected at the same time from the n element array. Now once you demodulate or remove the carrier effect from the n samples then the direction finding problem is same as frequency estimation problem mentioned in case of time domain data.

The above is true when there are less then 'n' frequencies present in the signals. Why this is true: the phase difference of the received signal on the array directly related to the angle of arrival. Kindly recall that phase of a sinusoidal is constant.  So time domain data with "fixed sampling frequency" is same as samples on Linear array with equal spacing element(fixed time) after demodulation .

Please understand it is true for narrow-band signal only not broad-band chirp like signal.

Hope it helps.


Effect of origin poles on stability? What will be stability if we have only one single pole at origin in s domain??
and what will be the case for multiple poles at origin in s domain?

A system with simple distinct poles on the imaginary axis (and note that the origin is on the imaginary axis) and no poles in the right half-plane is called marginally stable. If you have poles with multiplicity greater than $1$ on the imaginary axis, or if there are poles in the right half-plane, then the system is unstable.

For discrete-time systems, the same is true if you replace "imaginary axis" by "unit circle".

Transients in marginally stable systems do not decay, but neither do they grow without bounds. In practice, we usually want to avoid marginally stable systems. An ideal integrator is an example of a marginally stable system.


What features to use for multilayer perceptron model of an EEG signal? I have an EEG signal that I want to extract features to apply multilayer perceptron (MLP). What I should use, Fourier or wavelet coefficients?

Fourier analysis doesn't provide you with temporal information.  If you need temporal information you have several choices, two of which are the short term Fourier transform STFT, and wavelet transforms.

The STFT provides frequency coefficients in fixed time intervals across your data.  There is a trade off between the time and frequency resolution obtained.

Wavelet transforms pick apart features in your data at various time resolutions which may give you better time localization for your coefficients.  Additionally, you can choose different wavelet functions to control how features of your signal are extracted.  Another way of saying this is that your choice of wavelet function affects what your coefficients represent.

Most of the work on feature extraction that I've seen recently (this is in general, no specific to EEG signals) employ some form of wavelet analysis.

Here is an example paper: http://uprt.vscht.cz/vav/RSoftware/S7help.pdf


What does a star shaped autocorrelogram mean? Using the data in the upper graph, I get the weird autocorrelogram of the lower graph:-

I've never seen a autocorrelogram of this shape.  They're usually flat along the x axis.  Testing the code with simple sine waves produces the traditional flat bottomed graph. What causes the star shape, and is it valid or have I just done the processing badly? It seems specific to my data.

Python /numpy /scipy code:

correlation = signal.convolve(signalData, signalData[::-1], mode='full')

I'll go ahead and post a practical answer for you, and we can both see if someone can post a theoretical explanation.

The shape you see comes from the offset in your data.  If you were measuring voltage, I'd say you have an AC signal on a DC offset.

I don't know what you are measuring, but it spends most of its time at around 100.  It has short intervals where it changes drastically, and makes small changes almost constantly.

You have two ways to get rid of the triangle.

Take the average of all data points and subtract that average from all the data points.
Use a high pass filter to remove the offset.

If this were an audio signal, you could pick a lower cutoff frequency and be happy - there's some point where the audio frequencies are too low to be interesting.

Since you don't mention the source of the data (could be the sum of loans and borrows of elephants in an elephant lending library for all I know,) I'd go with the first option.

If it is some physical quantity, then you may want to consider the high pass filter, after all.

If the process generating the signal has some natural low frequency, then you could use that as your cutoff.

Or, you set the cutoff by the longest time period you want to consider.

If your data set represents 1 second, then autocorrelation for more than 1second is pointless - you could set the cutoff for one hertz.   At any rate, you would use 1/(time period) as the cutoff frequency.

If you have all the data in blocks at hand, you can use filtfilt to do your filtering.  That preserves the phase and the absolute timing of your signal.


Difference between delayed down-converted signal and a down-converted delayed signal I am having trouble understand what I think should be a pretty simple concept.  Conceptually I can understand that a signal that has been down-converted, and then delayed in time is different than a signal that has been delayed in time and then down-converted.

I am trying to compensate for this difference on an FPGA (I am down-converting the signal, delaying in time, and then up-converting)

The way I see it, the delayed (and attenuated signal) looks like:
$A e^{i(\omega_{\text{c}}-\omega_{\text{lo}})(t-\tau_1)} m(t-\tau_1)$

That is different from a down-converted delayed signal, so I need a phase correction constant:
$\phi = e^{i\omega_{\text{lo}}\tau_1}$

Which gives:
$A e^{i(\omega_{\text{c}}-\omega_{\text{lo}})t - i\omega_{\text{c}}\tau_1} m(t-\tau_1)$

So the correction is a function of the LO and the delay itself, right?  I have my sampled signal in I and Q, but I need to multiply it by the correction before up-converting, and that is where I am lost.

Let your original signal be (leaving out any constant multiplicative factor):

$$s(t)=e^{i\omega_ct}m(t)$$

Then the down-converted signal is

$$\tilde{s}(t)=e^{i(\omega_c-\omega_{lo})t}m(t)$$

and its delayed version is

$$\tilde{s}(t-\tau)=e^{i(\omega_c-\omega_{lo})(t-\tau)}m(t-\tau)$$

whereas the delayed and down-converted signal is

$$s(t-\tau)e^{-i\omega_{lo}t}=e^{i\omega_c(t-\tau)}e^{-i\omega_{lo}t}m(t-\tau)=
\tilde{s}(t-\tau)e^{-i\omega_{lo}\tau}$$

So the correction factor $e^{-i\omega_{lo}\tau}$ indeed depends on both $\omega_{lo}$ and the delay $\tau$.


How to de-noise raw sound data What techniques/algorithms can I use to remove noise from a raw recording of sound (voice)?

The purpose is to get a smoother graph (removing the "jaggedness"). What I have tried was to average-out small deviations by using the surrounding two values to either side, but ended up with distorted sound. What better way is there to "smooth-out" the graph?

For what you intend to do, a low-pass filter is the way to go. Your statement about filtering frequencies vs filtering amplitudes is incorrect. Your signal contains many components at many frequencies, the amplitude of which varies in time, and the high frequency components are those causing the "jaggedness" and you want to get rid of them. Not sure why you say your signal is "constant frequency" - maybe you are getting confused by the sample rate?

What you have tried (averaging) is indeed a special case of low-pass filtering, but one with a frequency response far from being ideal. You should try a properly designed IIR or FIR filter. In particular, the FIR filter is not very different from what you tried - this is just a weighted combination of the samples neighboring each sample. But the choice of the coefficients is important and ensures that only unwanted components are eliminated. Note that an FFT is not the way to go. This question crops up here quite frequently under different forms, but in short - FFT, messing with coefficients, IFFT - is a bad idea.

By design, the output of a moving average filter (what you implemented) has less energy than the input. It is thus impossible for a moving average filter to cause distortion. If the input signal is in the [-1, 1] range, there is no way for an averaging filter to yield values outside this range. The "distorted sound" you observed was probably due to an implementation error on your side (overflow / clipping of integer values, signed values treated as unsigned value, or maybe incorrect in-place processing)...

EDIT: one thing worth mentioning is that there are situations in which a speech signal actually has high frequency components (appears "jagged", is noisy) - for example during a sss or shhh ; and removing those with a low-pass filter will affect its brilliance. Ideally, you'll want your low-pass filter to be active only when you detect that your speech signal is voiced - and inhibit it when you detect an unvoiced, noisy consonant.