Does the inverse Fourier transform only produce stationary time signals? I have a question regarding the inverse Fourier transform and its relevance to non-stationary signals. And by non-stationary signal, I'm talking about a signals whose frequency content varies with time. I'll provide the example that made me originally ask this question.

Say for example I have a sinusoidal signal with a single frequency that goes on for ever. I could take the Fourier transform of this signal to obtain its frequency content. And because this signal is a single frequency that goes on forever, this signal is considered stationary.

However, now let's say I have a signal that is zero except for in some window in which the signal is sinusoidal. This nonzero window exists from t = 0 to 1 second. If I were to take the Fourier Transform of this signal over a finite duration (such that it acts on both the nonzero and zero part of the signal), I would get a frequency content that shows there exists a peak frequency content due to the nonzero window.

Now let's say I were to do the same thing, except this time the nonzero window exists from t = 6 to 7 seconds. If I were to take the Fourier transform, over the same finite duration, I believe I would yield the exact same frequency content (both phase and magnitude) from the earlier experiment.

Assuming what I said was correct, it should thus be IMPOSSIBLE to recover the exact signal from the Fourier transform since two different time domain signals have mapped to the same frequency domain signal. Am I correct with this statement? Furthermore, I believe the inverse Fourier transform will end up producing a stationary signal. As in, the frequency that was once contained just within the nonzero window will now be shown to be infinite from the inverse Fourier transform.

TLDR:

Does the inverse Fourier transform only produce stationary time signals?
Thus if I have non-stationary signal, take the Fourier transform and then take the inverse Fourier transform, I will now have a signal that is stationary?

The FT of a finite length sine wave burst consists of more than just the one frequency of that sine wave.  A finite width time domain signal has an infinite width FT.  Same for an FFT, a rectangular window produces a full width periodic Sinc in the FFT result.  So even if the phase of the 1 sine wave result bin stays constant, the phase of all the other frequencies will rotate as the position of the window changes.  Two different rectangular windows will NOT map to the same FFT, but convolve against any single sine wave differently.

Therefore an inverse FT or IFFT of the full spectrum can produce a signal whose statistics can vary significantly over different portions of the result, by reproducing a window which can be zero in parts of the result, and non zero in others.


Does image size impact compression efficiency differently in different formats? E.g., if an image is, say, 600px x 400px and the PNG format produces a smaller file size than JPG, does this hold true at any image size? If the image were 600,000px x 400,000px would the two file sizes still have the same relative proportions? Or is it possible that JPG's compression means that PNG will eventually produce a larger file size past some number of megapixels?

JPEG images can't be that big; their maximum size is (2¹⁶-1)×(2¹⁶-1). (doesn't really matter)
These are fundamentally different file formats, and their compression stems from different properties: Where JPEG is a lossy format who is able to drastically reduce the file size of content that has photographic properties (mostly: low-pass behaviour in 8×8 pixel blocks) at different levels of lossiness, PNG is lossless and compresses especially well on images with similar rows.
Thus, the size of the images depend on the content. You can't directly extrapolate from a small image to a larger version: it depends on the content of the original large image whether it compresses better in JPEG (for a given lossiness factor) or PNG.

Hence, no statement of the sort you seem to be looking for can be made, in general.

Specfically, for images that are large and that results in large swaths of identical pixels, JPEG always operates on blocks of 8×8 pixels. Even if there's a lot of monocoloured blocks, these still have to be stored and later entropy-coded in JPEG; this leads to repeating bit patterns, but not necessarily constant 0-bytes, having to go through a Huffmann-alike compressor.

In PNG, identical runs of pixels will be converted to differences to neighboring pixels; thus, large portions of 0-runs exist in large pictures with constant areas. The DEFLATE algorithm that follows this differential encoding is very efficient in reducing consecutive constant values to very few bits.


Real time numerical differentiation of signals I am trying to implement a PID controller in Python and I am having some problems with real time numerical differentiation of my my discrete signal.

I am using the following method:

$$d{\frac{x_n}{t}} = \frac{x_{n-1}-x_n}{d_t}$$

Where $d_t$ is the time difference of the measurement time of $x_n$ and $x_{n-1}$.

The differentiation is accurate to some point but the results on real time systems are not as accurate as I would like to be. The differential controller does not really improve the stability of the system in most cases. At high gains it starts introducing noise into the system.

Most numerical differentiation methods recommend going into the future ($x _{n+1}$) and sadly for my system that's not possible. Such as:

$$\frac{x_{n+1}-x_{n-1}}{2\times d_t}$$

Are there any methods for calculating real-time differentiation with better accuracy?

Your original differentiator, which should be $x(n)-x(n-1)$, is called a "first difference" differentiator. That differentiator amplifies high-frequency noise. As a next step I suggest you try what's called the "central difference" differentiator defined by:

$$ \mathit{Diff} = \frac{x(n)-x(n-2)}{2} $$

which does not amplify high-frequency noise.


recover image by only magnitude of image fourier transform I have a image and I want calculate fft2 of it, after it I want recover image only by magnitude of it. how can I achieve this work?

I=imread('1.bmp');
I = rgb2gray(I);
imafft = fft2(double(I));

 % Gets magnitude and phase.
 mag = abs(imafft);

 % Generates the modified spectrum
 newfft = mag.*exp(0*angle(imafft));
 ima_out = ifft2(newfft);
 imshow(abs(ima_out));

this is return me white page.

Check out this classic example from Oppenheim, A. V., & Lim, J. S. (1981). "The importance of phase in signals". a) and b) are the original images, c) is the image created using the phase of a) with the magnitude of b), d) is the image created using the phase of b) and the magnitude of a). Phase carries most of the information in an image.


What is different between Euclidean distance and the Geodesic distance? We have two kinds of distance in image: Euclidean distance and the Geodesic distance. What is different between them? Could you show mathematic formula or visualization to make it clear? The reference said

    Difference between the Euclidean distance and the Geodesic distance calculated with the Fast Marching Method. The Geodesic distance is the distance of the minimum length inside the figure path and the Euclidean distance is the straight line distance

This is example

Simple: the Euclidean distance completely ignores the shape when finding a path from the start point to the end point while, for the geodesic distance, the path is constrained to be within the given shape.

That's why the distances at the bottom left of the figure are so different.


Upsampling in Python I'm trying to perform upsampling in Python. Given the signal 's' and factor of upsampling N, I wrote the following code:

> y = np.zeros(N*len(s))
> for i in range(0, N*len(s), N):
>     for j in range(len(s)):
>             y[i] = s[j]
> return y

However, when I have a signal with first value '1' and factor of upsamling 3, my output of the function is like [1, 0, 0, 1, 0, 0] etc, no matter what the rest of the signal looks like. Is there a way to make this code work properly?

I think you are trying to do this:

y = np.zeros(N*len(s))
for i in range(0, len(s)):
    y[i*N] = s[i]

return y

There are other ways in python that avoid explicit loops, for example:

y = [0] * (N * len(s))
y[::N] = s


Low recognition rates using GMM for image classification I want to use GMM for image classification.
So, I have extracted SIFT features from each image in the corpus.
Then, I apply EM algorithm to learn GMM parameters (I have coded it in matlab).
I get very low Recognition rates for unseen images when I test the model on larger number of classes (>=4).
Who can explain this result ?  Thanks.

Are you trying to model the distribution of the SIFT descriptors using GMM? Are you using the actual 128-element descriptors?  If so, then you are trying to model a distribution in a 128-dimensional space, which is not likely to work well.

You need to reduce the dimensionality. If you insist on using GMM, then the simplest thing to try is to use principal components analysis (PCA) to reduce the dimensionality of your descriptors.  I would not use more than 10 dimensions.  Here is a paper, that tried kernel density estimation on reduced SIFT descriptors among other things.

Alternatively, you could try other algorithms, such as Bag Of Features, or the Pyramid Match Kernel.


Normalized cross-correlation in frequency domain I never worked with signal processing and never really used Fourier transforms before, still I am working on a project consisting on taking the output of an accelerometer to detect some movement features between touching surfaces. Using a frame of the output, I compute the FFT. After analyzing the data for multiple experiments, it appeared to make sense to use gaussians to isolate bands of frequencies and use that to classify the features.

That does not work well, however, if I change the sensor I am using, the sampling rate or frame size. Thus, I want to calculate the normalized cross-correlation. Problem is: I could not find an algorithm to calculate the NCC in the frequency domain and I did not figure out how to transform the translated gaussians to the time domain.

For instance, I have a frequency spectrum from -F/2 to F/2 and L samples. I converted that to L/2+1 frequency bins from 0 to F/2 and defined a band using a gaussian with standard deviation s and centered on c.
To get the cross-correlation between the spectrum and the gaussian, I simply do the sum of the element-wise multiplication of them. Now, if I want that normalized, what can I do? Is there a way to normalize it right from the frequency domain?

I was thinking on the math behind it and concluded that calculating the NCC in the frequency domain is basically the same process as in the time domain:

$F[k]$ is the (discrete) Fourier transform of my sample, $G[k]$ is the gaussian curve in the frequency domain. The cross-correlation is
$$CC = \sum^{L-1}_{k=0}F[k]{\cdot}G[k]$$
Then, to get the normalized cross-correlation we first normalize the vectors $F$ and $G$, as would be done to any vector:
$$\hat{v} = \frac{v}{|v|} = \frac{v}{\sqrt{\sum^{L-1}_{k=0}v^2[k]}}$$
Thus, the normalized cross-correlation should be the dot-product of the normalized vectors $\hat{F}$ and $\hat{G}$:
$$NCC = \frac{\sum^{L-1}_{k=0}F[k]{\cdot}G[k]}{\sqrt{\sum^{L-1}_{k=0}F^2[k]\cdot\sum^{L-1}_{k=0}G^2[k]}}$$


Compensating effects of a system with a known transfer function Suppose we have a system which we want to know the exact transient times. In ideal case, we can extract the transient times, but in practice it will be affected by another system with a known transfer function. We know  the system with a transfer function  changes the transient times and also the shape of output. How can we obviate its effects on the (main) system?

That is pretty broad question, as it depends on the system.

However, your problem is extremely common: Every communication system works because a receiver is able to estimate the original events (symbols sent by the transmitter) after the signal has gone through a system with transfer function (the channel).

So, receivers need to do exactly that: reverse the effect of the channel as good as possible. The thing that does that is usually called an equalizer; the equalizer you pick, and how much info about the system you need to make that equalizer work, depends on your "main system".

So, if you need to solve that: go ahead; define your main system as well as possible. Find stochastic descriptions for the parts of the system that you don't know with certainty. Then, find a channel model that fits that model, and look for equalizers used with such a channel.


Relating phase shift to the ROL/ROR operations Is phase-shifting in continuous-time analogous to the rotate-left (ROL) and rotate-right (ROR) operations in discrete-time?

If you mean phase shift of a sinusoid, then yes, a phase shift is equivalent to a translation.


SNR in one-dimensional additive channel with uniform noise and binary antipodal input A one-dimensional addictive noise channel, $Y=X+N$, has uniform noise distribution $N \sim U(-\frac{L}{2}, \frac{L}{2})$ and binary antipodal input constellation with equally likely input values $x=\pm 1$. Find $\rm SNR$ as a function of $L$.

I know that the $\rm SNR = \frac{\rm Signal \, \,Power}{Noise \, \,Power}$, yet I couldn’t figure out how to find these power. Any hints?

Signal power $=E[X^2]$; Noise power $=E[N^2]$. I assume that you know how to calculate these expectations.


Calculate displacement from acceleration signal data using Fourier transformation I recorded acceleration data from an accelerometer attached to a vehicle, and I want to calculate displacement using the Fourier Transformation Integration method.
I used software called vibration Data ToolBox.  This software has a built-in Fourier transformation which lets me save transformed acceleration data into a complex number (real + imaginary).
According to the rules of the Fourier transform, I can perform integration on the time-domain data by dividing the transformed acceleration data by the scale factor (-2πf) where f is the frequency in Hz.
However, the transformed data is a complex number.  Should I divide the real part or the imaginary part or the magnitude by the scale factor?
Also, the software can perform Inverse Fourier transformation just for complex transformed data to obtain displacement in the time domain. how to calculate displacement?
I attached transformed data as Excel sheet.Excel sheet

I can perform integration on the time-domain data by dividing the transformed acceleration data by the scale factor (-2πf) where f is the frequency in Hz.

Well, not quite, but it looks like you have pieces of it.  To integrate once, you divide by $2i \pi f$ -- note the fact that the divisor is purely imaginary.  To integrate twice you just divide by that factor twice.  You can do that explicitly, or you can note that $\left (2i \pi f\right)^2 = -4 \pi^2 f^2$

However, the transformed data is a complex number. Should I divide the real part or the imaginary part or the magnitude by the scale factor?

The key phrase here is "complex number".  Almost always with complex numbers, you do arithmetic on the whole thing, real and complex parts equally.  If you're just plucking out the real or the imaginary parts out to do math on, you're almost certainly taking some shortcut, and you should understand the "right" way to do that task.
In this case, you're just doing plain old arithmetic, so you should divide the entire number (meaning: the real and the imaginary parts) by $-4 \pi^2 f^2$.

Also, the software can perform Inverse Fourier transformation just for complex transformed data to obtain displacement in the time domain. how to calculate displacement?

Theoretically, divide throughout by $-4 \pi^2 f^2$ and do an inverse Fourier transform.  In reality, you'll have some non-zero signal content at low frequencies, and dividing that by 0 will cause problems.
So you need to multiply it by some function $A\left(f\right)$ that is $A\left(f\right) = 0$ for $f = 0$, $A\left(f\right) = -\frac{1}{4 \pi^2 f^2}$ for $f$ above some threshold frequency (which will be determined by your accelerometer, BTW), and that transitions smoothly in between.
At worst, just make a function $$A(f) =
\begin{cases}
0 && |f| < f_0 \\
-\frac{1}{4 \pi^2 f^2} && |f| \ge f_0
\end{cases}$$
This will hopefully work well enough to tell you that you're doing your work right, but will suffer from time-domain ringing around $f_0$.  Then you can either try out ways to make a smooth transition, or ask back here and we'll suggest something.


Impulse response from difference equation without partial fractions I have a system with input $x[n]$ and output $y[n]$, described by the difference equation

$$y[n] - 0.3y[n-1] + 0.1 y[n-2] = 3x[n] - x[n -1]$$

and I am trying to find the frequency and impulse responses. The frequency response I found was

$$H(e^{jw}) = \frac{3-e^{-jw}}{1 - 0.3e^{-jw} + 0.1e^{-j2w}}$$

but I am having trouble  converting this to $h[n]$ because the denominator can't be factored using partial fractions.

Also, I tried using the recursion method to find $h[n]$ but I am unsure whether the system is causal or not, so I don't know where to start for the recursion. Could anyone give me some tips on what to do next?

First of all you can apply partial fraction expansion method to get the inverse DTFT to find the impulse resposne h[n] for this system. You have to assume casual ot anti-causal system.

$$H(e^{jw}) = \frac{3-e^{-jw}}{1 - 0.3e^{-jw} + 0.1e^{-j2w}} = \frac{3-e^{-jw}}{ (1 - \alpha e^{-jw})(1 - \alpha^* e^{-jw})}$$

$$H(e^{jw})= \frac{A}{1 - \alpha e^{-jw}} + \frac{B}{1 - \alpha^* e^{-jw}}$$

Where $\alpha = 0.1500 + j0.2784$ for this particular case. Hence assuming causality and using inverse DTFT look-up tables, deduce as:

$$\boxed{ h[n] = A \alpha^n u[n] + B (\alpha^*)^n u[n] }$$

Time domain recursion is also possible, requiring only two auxiliary values.
For example assuming causality one can use recursion from the LCCDE to obtain $y[n]$ for all $n \geq 0$ from the given LCCDE
$$y[n] - 0.3y[n-1] + 0.1 y[n-2] = 3x[n] - x[n -1]$$

by first rearranging it for $y[n]$ as
$$y[n] = 0.3y[n-1] - 0.1 y[n-2] + 3x[n] - x[n -1]$$

Then the recursion follows for $n \geq 0$ assuming that $x[n]=\delta[n]$
$$
\begin{align}
y[0] &= 0.3y[-1] - 0.1 y[-2] + 3x[0] - x[-1] \\
y[1] &= 0.3y[0] - 0.1 y[-1] + 3x[1] - x[0]  \\
y[2] &= 0.3y[1] - 0.1 y[0] + 3x[2] - x[1] \\
y[n] &= 0.3y[n-1] - 0.1 y[n-2] \\
\end{align}
$$

Note that in computing $y[0]$ you need $y[-1]$ and $y[-2]$ as auxiliary (initial) conditions. Also note that $x[-1]=0$ and $x[0]=1$ and $x[n]=0$ for all remaining $n>0$.

A similar analysis fallows for the anti-causal case.


Creating an array whose DFT is all whole numbers I'm working on a simple FFT implementation and currently running some unit tests. Long floating point numbers in C and assembly become very hard to keep track of, and this whole thing is error prone. I was able to some up (purely experimentally) with a few 4-tap and 8-tap arrays whose DFT is all whole numbers. I'm only testing DFT arrays of length $2^N$, and I could, of course, repeat my shorter signals several times until I achieve desired length, but this would give me a whole bunch of zeros, and that's not a good thing to test for. In summary, I need to come up with signals that satisfy the following:

Signals are of length $2^N$
Coefficients of both the signal and its DFT must be whole numbers
Signals are real
Signals must not consist mainly of zeros

Could there be a general way of coming up with such things?

Update
To explain number 4 above, let me demonstrate what I mean with an example. I don't want a situation in that was produced by repeating something that has a known integer DFT several times.DFT of [1 5 -1 0] is [5 2-5i -5 2+5i]. If you simply repeat [1 5 -1 0] four times to get

[1 5 -1 0 1 5 -1 0 1 5 -1 0 1 5 -1 0],

its DFT would be

[20 0 0 0 8-20i 0 0 0 -20 0 0 0 8+20i 0 0 0].

This shows that repeating something many times does not add any new non-zero integer numbers to the DFT. This procedure simply inserts zeros between existing DFT coefficients (and scales them depending on definition of the transform). This isn't of much help. So zeros are welcome, but not the zeros achieved this way.

The DFT matrix has $e^\frac{-2 \pi j}{2^N}$ as coefficients. For $N >= 3$, these coefficients might have irrational numbers as real/imaginary parts (eg $\frac{\sqrt 2}{2}$ for $N=3$); so all these irrational terms will have to cancel out exactly in the DFT computation. There are ways of achieving that (repeating is one of them) ; but this keeps you away from the general case that your code has to test! For example, this will create symmetries in your data that might make an error in the bit reversal, or in the higher orders (3rd level or above) of butterflies undetected!

Another thing to take into account is that due to the pure floating point implementation, even if the expected result should be a whole number, your FFT code, especially on large sizes, might find a result like 0.99999999 (Heck, some whole numbers are not even representable as floats!) ; and this wouldn't be incorrect. Plausible on doubles, very likely of floats.

So you'd rather use very generic data (ie: that will yield some irrational numbers) ; and some unit tests that guarantees you that your results are within some tolerance of values computed by a reference implementation such as FFTW. For example in gtest there is EXPECT_FLOAT_EQ or EXPECT_NEAR to check for float equality (up to representation errors) and for float equality given a tolerance.


Gain at critical frequency for Butterworth and Chebishev filters In Chebyshev filter the critical freq occurs with gain Gp But in butterworth filter it is always at -3dB; why does this happen? Can anyone explain intuitively on this?

If I understood your question correctly, Butterworth filter is maximally flat (no ripples) in the passband frequencies, and at the cutoff frequency, the gain is down by 3dB (3.01 to be exact). But disadvantage is slow roll off from passband. See https://en.wikipedia.org/wiki/Butterworth_filter especially the figure showing magnitude roll off.

For Chebyshev filter type 1 filter, the roll off is more sharp but result is ripples in pass band. You can control this ripple, which will in turn vary the Gp that you mentioned. See https://en.wikipedia.org/wiki/Chebyshev_filter. If you specify $\delta$ as the ripple factor in dB scale (like 1dB), $\epsilon = \sqrt{10^{\delta /10}-1}$ in linear scale. The gain Gp = $\frac{1}{\sqrt{1 + \epsilon^2}}$ in linear scale. The sharpness of your passband and passband ripple specification will affect $\epsilon$ which will inturn affect Gp.


Comparing multiple signals for similarity I have multiple (between 2 and 100) signals and need to determine when a significant number diverge from the rest. We're exploring machine learning techniques, but we also want tackle this as a signal processing problem and see where we get the best results.

This very informative post suggests that best results come from a weighted ensemble of techniques, including:

Similarity in time domain (static): Multiply in place and sum.
Similarity in time domain (with shift*): Take FFT of each signal, multiply, and IFFT. (matlab's xcorr)
Similarity in frequency domain (static**): Take FFT of each signal, multiply, and sum.
Similarity in frequency domain (with shift*): Multiply the two signals and take FFT. This will show if the signals share similar spectral shapes.
Similarity in energy (or power if different lengths)

But this is a fairly high-level outline. Can anyone point me to a more thorough discussion of these techniques, preferably with some python code or in lieu of that, some code in R?

...best results come from a weighted ensemble of techniques...

Maybe they do, depending on the application. But each one of the similarities mentioned, is equivalent to the other at least when we are talking about signals originating from linear systems.

Cross correlation provides very good results especially if you are trying to figure out if a signal is noise or not and a very good application of that is in the rake receiver. And via the correlation's link to the Discrete Fourier Transform (DFT) that enables us to carry it out in the frequency domain we can start to see the equivalences between the given statements.

A more informative division would be between linear "metrics" and non-linear metrics. But, the field of "similarity metrics" is vast, so, it is best to start with what is known from the problem and then work towards a solution.

What I am going to do here is to talk a little bit about some metrics of similarity which you might find useful but then also write a little bit about other ways of doing the same thing that don't require the assessment of similarity as I think that that second part is closer to what you are trying to achieve (but not what you are asking directly).

On Similarity Metrics:

Cross-correlation will work fine if your signals have been consistently sampled to each other. That is, each device's sampling frequency remains absolutely stable but the relative points in time that the recording started can be varied. This will insert a delay between waveforms which is very easy to recover by doing multiple cross-correlations at different delays.

If you want to look at specific frequencies, you can either filter the signals prior to cross-correlation or use coherence.

Another thing you can try is Principal Component Analysis. To do this, pack all the signals that you receive from $N$ devices into a matrix and send it to PCA. PCA will try to decompose these signals into a smaller set of "principal" (uncorrelated) signals and will return a matrix with those elementary signals and their strengths. It does this using covariance, which is very closely related to correlation. The point here is that if all of the signals are the same, there are no "components" and the vector of strengths is going to be zero. But, if there is significant variation within the signals, the vector of strengths will start exhibiting some structure. You can then use the vector of strengths in a simple classifier that discriminates whether the indicated variation is negligible or worthy of an alarm.

Now, it might be that the sampling period of each GPS device is not stable (very likely). And this now means that if you were to put "one signal on top of another" (as if they were printed on transparency slides), it would be possible to align some peaks, but not others.

In that case, you could use Dynamic Time Warping (DTW). Imagine that you listen to a vinyl record with a lot of wow. You can remove that wow, by speeding up the recording when the record player slows down and vice versa. This is exactly what DTW does. You pass to it two signals and it returns a matrix that  tells you where to slow down and where to speed up. After you apply this transformation to one of the signals, the peaks would align perfectly. The total amount of deviation in this matrix is a metric of similarity. But, DTW is not exactly fast...particularly if you have to consider $\frac{N*(N-1)}{2}$ unique combinations (because all of these metrics are symmetric) between all 2 out of $N$ signals.

On the specific application:

So, the application is a GPS interference alert. There are two use-case scenarios that I can see at the moment:

GPS Jamming

Where the receiver cannot access the signals from the satelites because a local source is "flooding" its input with "noise".

GPS Spoofing

Where the receiver thinks they are listening to a satellite, but actually they are listening to a loud black box that behaves like a satelite.

Since it is impossible to raise an interference alarm based solely on information available to a device, the idea is to incorporate information from other data sources (nearby devices / data) in order to confirm the "story" that the GPS receiver on this device is listening to.

How can we get a gauge on the quality of the fix on a given device?

There are two obvious choices:

The Dilution of Precision
The GPS Quality Indicator Field.

Between these two, the Dilution of Precision is a "softer" indicator (wider range of values) of the quality of the fix.

So, assuming that there is an infrastructure for collecting GPS data from a set of devices (like Waze or Glypmse do for example), then the following alarms can be raised:

To counteract GPS Jamming:

Given the timestamp and the GPS Quality Indicator for a set of devices

Select a point in time
Select the subset of devices which have provided "bad" GPS quality data at or around a reasonable time interval around the selected point in time.
Cluster the provided GPS locations in space

This can be done ridiculously fast
A Threshold is required to define which two clusters are considered apart enough. This can be discovered automatically or be set to a reasonable value. Bunches of receivers approximately 1km apart should be considered separate, for example.

Assess the "size" of the "patch" of bad reception.

An obvious choice here would be to take the centroid of the cluster and find the furthest data point (so, one radius), or find the mean radius of the points or take their convex hull or use another more accurate algorithm that describes the points.

Decide if the patch is big enough to warrant raising the alarm.

Obviously, lots of receivers in a short tunnel or underpass or close to tall buildings are not supposed to raise alarms.

To counteract GPS Spoofing:

This is much more challenging because the problem there is to identify an unjustifiable deviation. This assumes that you have some way of knowing what is "correct" and that is the key difficulty. The easiest thing to tackle here is the GPS being spoofed over a wide area (more than one device).
In this case, you want to identify a cluster of devices that start behaving consistently all together. As you can understand, if 10 sailing autopilots were listening to the same spoofed signal and we start rotating the frame of reference then all of them will be taking consistent decisions at the same time. The above algorithm more or less remains the same but this time you want to cluster devices whose bearing (derived bearing) deviated by a consistent value at the same time. That is, they all appear to be adjusting at the same time and by the same amount. (As you can see here, no reason to assess the similarity of specific signals, we work with location directly).
If you are trying to understand if a single device's GPS is being spoofed, then that is even more challenging and there you would need to have access to other data as well (again to try to confirm the story the GPS is telling you).

Hope this helps


Faster R-CNN and 4-step Alternating Training I am reading the paper Faster R-CNN available here. In section 3.2, page 6 bottom left, the authors describe what they refer to as a 4-step Alternating Training between the Region Proposal Network, used to obtain region proposals used in the Fast R-CNN method for object detection. They write:

  In this paper, we adopt
  a pragmatic 4-step training algorithm to learn shared
  features via alternating optimization. In the first step,
  we train the RPN as described in Section 3.1.3. This
  network is initialized with an ImageNet-pre-trained
  model and fine-tuned end-to-end for the region proposal
  task. In the second step, we train a separate
  detection network by Fast R-CNN using the proposals
  generated by the step-1 RPN. This detection network
  is also initialized by the ImageNet-pre-trained
  model. At this point the two networks do not share
  convolutional layers. In the third step, we use the
  detector network to initialize RPN training, but we
  fix the shared convolutional layers and only fine-tune
  the layers unique to RPN. Now the two networks
  share convolutional layers. Finally, keeping the shared
  convolutional layers fixed, we fine-tune the unique
  layers of Fast R-CNN. As such, both networks share
  the same convolutional layers and form a unified
  network.

What I do not understand is the third step. They say they use the detector network to initalize RPN training. I am not sure what this means. The detector network outputs bounding boxes and class labels for objects within a given image.

Consider the layers specific to RPN. We have a certain number of convolutional layers, followed by sliding windows with anchor boxes being fed to regression and classification (class label layers). Are the authors saying that the initially trained convolutional layers are now discarded for RPN and replaced with those trained by Fast R-CNN in step 2? Or do they mean something else.

Since no one answered this, I will post my attempt after further research.
Here is my understanding after having thought about it more. Lets say we are using VGG-16 as the backbone.
There are three fundamental components, the backbone, the RPN heads (objectness and RPN regression), and the ROI heads (classification and regression). Using the RPN heads, we can sample positive and negative ROIs.
Step 1) Train a region proposal network (backbone + RPN head) with the VGG-16 backbone pre-trained on imagenet with the head popped off. Some layers of the backbone will be trainable.
Step 2) Using the object proposals (300 or so per image) generated from the network obtained in step 1) (the backbone trained with the rpn), train a second network with another separate VGG-16 backbone pre-trained on imagenet, and a network head as in fast R-CNN (we will call ROI head). Then as in fast R-CNN we train for object classification of a region, as well as regression to learn the displacement from the proposal to a bounding box around the object (if the highest activation in the softmax does not correspond to a "background" class, i.e. no object).
Step 3) Use the backbone from step 2) for RPN training, which will be fixed (= not trainable), with the RPN head trained from step 1), and fine tune the RPN head.
Step 4) Finally keeping the backbone from step 2) fixed, and the RPN heads from step 3) fixed, fine tune the ROI Heads from step 2) to train bounding box regression and classification once more on the region proposals generated from
the backbone and RPN head.
Edit: Original answer had many errors so I came back to correct them, also original question doesn't really make sense now that I read it (I was new to this topic). Please just read this answer.


White gaussian noise analysis deduction I´m stuck in a deduction analysis of the variance of a gaussian white noise signal in a "integrate-and-dump detector" of a baseband data transmission receiver, where $n(t)$ is white noise with double-sided power spectral density $N_0/2$ [W/Hz]

I can understand all the steps except when they deduce

How do you get to this last expression?
Thank you.

They say that $n(t)$ is white noise with a double-sided power spectral density (PSD) of $N_0/2$, i.e., the PSD is given by
$$S_n(f)=\frac{N_0}{2}\tag{1}$$
The auto-correlation function is the inverse Fourier transform of $(1)$, which is
$$R_n(\tau)=E\big\{n(t)n(t+\tau)\big\}=\mathcal{F}^{-1}\big\{S_n(f)\big\}=\frac{N_0}{2}\delta(\tau)\tag{2}$$


Demonstrate the mean of the sample variance Let's suposse I have a vector of elements $x(n) = \{x(0), x(2), \cdots ,x(N-1)\}$ from a random process X of mean $\mu_x$ and variance $\sigma_x^2$. I want to see if I can estimate the mean and variance of the random process with my finite vector.

So, we say an estimator is unbiased if the mean of the estimator is the same as the mean of the random process. For example, the sample average is:

$$\hat{\mu}_x = \dfrac{1}{N} \sum_{n = 0}^{N-1} x(n)$$

Taking average in each member:

$$E[\hat{\mu}_x] = \dfrac{1}{N} \sum_{n = 0}^{N-1} E[x(n)] = \dfrac{1}{N} \sum_{n = 0}^{N-1} \mu_x = \mu_x$$

Thus, the sample mean is unbiased.
My problem comes with the mean of the sample variance. The definition of the sample variance is the following:

$$\widehat{\sigma_x^2} = \dfrac{1}{N} \sum_{n = 0}^{N-1} (x(n) - \hat{\mu}_x)^2 $$

I have to prove that:

$$E\left[\widehat{\sigma_x^2}\right] = \sigma_x^2 \left(1 - \dfrac{1}{N} \right) $$

I have tried to expand the square parenthesis:

$$\widehat{\sigma_x^2} = \dfrac{1}{N} \sum_{n = 0}^{N-1} E[x^2(n)] + E[\hat{\mu}_x^2] - 2 E[x(n) \hat{\mu}_x]$$

Nevertheless, when I have to calculate $E[\hat{\mu}_x^2]$ (I suppose it's this way):

$$E[\hat{\mu}_x^2] = \dfrac{1}{N^2} E\left[ \left( \sum_{n = 0}^{N-1}x(n) \right)^2 \right]$$

I don't know how to do it. Moreover, I don't know if $E[x(n) \hat{\mu}_x] = E[x(n)] E[\hat{\mu}_x]$ (I think that's not correct in general), and I have tried to solve it in another way, but anything comes to my mind.

I hope someone can help me. Thank you for your responses.

The key to prove
$$
E\left[\widehat{\sigma_x^2}\right] = \sigma_x^2 \left(1 - \dfrac{1}{N} \right)
$$
is to add and subtract $\mu_x$ in the expression for the sample variance before working out the square. A detailed proof is available on the Wikipedia page for Bessel's correction.


If I line in the sound from a video game console, can I detect when a sound is played? The question is as the question is. It won't be this isolated sound though, it'll be a sound on top of background music. I think I can get a hold of the sound file itself from the data, but at worst I'd just make my own recording and ask for a less accurate match - if variation is possible in whatever method may be proposed.

Theoretically you could separate your desired sound by ICA method or other blind source separation methods from the background music, then design a matched filter by your separated sound. The matched filter have high response when the desired pattern (which is your desired sound) exist in the input even when it buried in noise (background music).

The matched filter has acceptable performance but i'm not sure about the ICA and other source separation methods and highly depend on your input signals (desired and background sound).


Quadrature Amplitude Modulation: I and Q Same I am trying to understand a previous design that uses Amplitude Modulation.  Instead of just using the real part of the signal, the designers put the voice signal on both I and Q.  Both I and Q have the same data.

How would this affect the spectrum of the resulting AM signal and how would this affect the performance of the demodulator (assuming we are demodulating using a simple envelope detector)?

It doesn't really affect anything at all. If you modulate a signal on the I channel only, then it would look something like:

$$
x_1(t) = m(t) \cos(2\pi f_c t)
$$

If it were on the Q channel only, then it would look like:

$$
x_2(t) = m(t) \sin(2\pi f_c t)
$$

$x_1(t)$ and $x_2(t)$ only differ by a phase shift of 90 degrees on the carrier. Likewise, if you were to split the signal evenly in power across the I and Q channels, it would look like:

$$
x_3(t) = \frac{1}{\sqrt{2}} m(t) \cos(2\pi f_c t) + \frac{1}{\sqrt{2}} m(t) \sin(2\pi f_c t)
$$

Again, in the complex baseband form that is often used for analyzing signals, this corresponds to just a phase shift of 45 degrees from the original signal $x_1(t)$. In most systems, you will have a random phase shift between the transmitter and receiver anyway; this should not have any effect on the performance of your receiver.


How to evaluate CCDE at a given frequency? Consider the following Constant Coefficient Difference Equation (CCDE) given below

$$y[n] - \frac{1}{2}y[n-1] = 2x[n] - 5x[n-1] - x[n-2]$$

Let $H(e^{j\omega})$ denote the transfer function of this system. What is $H(e^{j\pi})$?

You need to follow the following steps (reponses at each steps are in a spoiler block, try to find them by yourself before checking). I will use the usual notation $z=e^{j w}$ for simplicity.

Find the relation between the transfer function $H(z)$ and the discrete-time Fourier Transform of the input $X(z)$ and the output $Y(z)$. This should be a basic results of your course

   $H(z) = \dfrac{Y(z)}{X(z)}$

Take the (discrete-time) Fourier Transform of your CCDE. For that, you will need to use two basic properties of the TF: linearity and time-shifting

   $Y(z) - \frac{1}{2} z^{-1} Y(z) = 2 X(z) - 5 z^{-1} X(z) - z^{-2} X(z)$

Compute $H(z)$ from the previous equation

   $H(z) = \dfrac{2 - 5z^{-1} - z^{-2}}{1-\frac{1}{2} z^{-1}}$

Finally take the value for $z = e^{j \pi} = -1$

   $H(e^{j \pi}) = \dfrac{2 + 5 - 1}{1+\frac{1}{2}} = \dfrac{6}{\frac{3}{2}} = 4$

This last result mean that your filter, at Nyquist frequency, has

   no phase delay and a gain of 4.


Understanding the transfer function of an FIR filter I'm currently studying FIR filter and am having trouble understand how the following equation works, and it's implication.

$$ y[n] = h[n] * z^n = H(z) \cdot z^n $$

I don't really understand how this works given that

$$ y[n] = h[n] * x[n] $$
$$ Y(z) = H(z) \cdot X(z) $$

It seems to me that $ x[n] = X(z) = z^n $.  What's going on here?  What domain is $ z^n $ in?

I thought that convolution is done on two functions in the n-domain; however, in the equation $ y[n] = h[n] * z^n $, $ z^n $ does not appear to me to be in the n-domain.

If it is in the n-domain, then I don't understand how the other part works, $ y[n] =  H(z) \cdot z^n $, because I thought that then it should be

$$ y[n] = Z^{-1}(\,H(z) \cdot Z(z^n)\,) $$

Thanks.

Edit:
Here's the page from the book with the equation.  It does seem to me that $z^n$ is used in both the time and frequency domain.

Looks like a notation or transcription error. The "z" that's used as the transform variable to calculate the Z-Transform can not show up as an independent variable in the time domain as well. One of them needs a different name.


Is a "Spectral Tilt Filter" the best way to get an arbitrary slope and cutoff LPF? I need an arbitrary slope audio LPF for which both the cutoff frequency and slope in dB/oct can be smoothly modulated in real time.

Is the "Spectral Tilt Filter" described by Julius Smith the correct or best way to do this?

https://ccrma.stanford.edu/~jos/spectilt/

https://ccrma.stanford.edu/~jos/spectilt/spectilt.pdf

If so has anyone implemented this in C++ before or would it be hard to do?

yeah, that's the way we've been doing pink noise filters (or "pinking filters") since the 80s.

if you ignore edge effects it's made by alternating these pole functions:

$$ \log|H_{2k}(j\omega)| = -\tfrac12 \log\left( 1 + (\tfrac{\omega}{p_k})^2 \right)  $$

and these zero functions:

$$ \log|H_{2k+1}(j\omega)| = \tfrac12 \log\left( 1 + (\tfrac{\omega}{q_k})^2 \right)  $$

and the total frequency response is

$$ \log|H(j\omega)| = \sum_n \log|H_{n}(j\omega)| $$

(this is not done yet)


Cross power spectral matrix from PSDs and coherence using matrix notation I have a vector containing all the PSD of a series of signals

$\underline{s} = \left[ |S_{x_1x_1}| |S_{x_2x_2}| |S_{x_3x_3}| ...\right]^T$

I also have a matrix $\gamma^2$ containing all the coherences

$\gamma_{ij} = \frac{|S_{x_ix_j}|^2}{S_{x_ix_i} S_{x_jx_j}}$

I would like to obtain the matrix of the cross-spectral densities, i.e. a matrix $\mathbf{S}$ so that
 $S_{ij} = |S_{x_ix_j}|^2$

Is it possibile to obtain it from $\underline{s}$ and $\mathbf{\gamma}$ using matrix notation?

The only way I found at the moment is:

$\mathbf{S} = (s s^T) \circ \gamma$

where $\circ$ indicates the Hadamard (element-wise) power and multiplication. This is however seldom used and usually element-wise multiplication is avoided in linear algebra.

Is there another way to do it?

$S= diag(\underline{s}) \cdot \underline{\underline{\gamma}} \cdot diag(\underline{s})$


How do I deconstruct a signal with fast fourier transformation in Python? I want to use the fast Fourier transformation in python to deconstruct a periodic signal into its harmonic periodicities.

Since the FFT of the original signal is the same as adding the FFT of the periodicities, I was calculating the power spectrum of the FFT of the original signal and then tried to extract the several signals. You can see my code for the extraction below:

def harmonic_period_calc(old_fft, peak_index, peak_value):
"""
    description: this function recalculates the fft so that only one peak in
                the psd should be visible

                 furthermore it uses the ifft to calculate a harmonic signal,
                 from which you are able to calculate the time indexes via
                 its peaks

    old_fft: the old fourier tansformation of our signal
    peak_index: the index, of peaks in the power spectrum
    peak_value: scalar-factor, with which we are going to shorten the other
            values with

"""

    new_fft_one = old_fft[:int(len(old_fft)/2)]
    new_fft_two = old_fft[int(len(old_fft)/2):]

    new_fft_one[np.delete(np.arange(new_fft_one.shape[0]), peak_index)] *= 1/peak_value
    new_fft_two[np.delete(np.arange(new_fft_two.shape[0]), -peak_index)] *= 1/peak_value

    new_fft = np.concatenate((new_fft_one, new_fft_two))

    period = scipy.fftpack.ifft(new_fft)

    return period

Unfortunately this is only working okay. Because the peaks of the inverse FFT is slightly off from the original period.

So my question is, does anybody know, how I will be able to do this?

If necessary, this is how I calculate the FFT and the power spectrum:

temp_fft = scipy.fftpack.fft(signal)
temp_psd = np.abs(temp_fft) ** 2
fftfreq = scipy.fftpack.fftfreq(len(temp_psd), 1/len(signal))
i = fftfreq > 0

Also the technical information from me:

Python 3.7
scipy 1.3.1
pandas 0.25.3
numpy 1.17.2

If you have a periodic signal, that is a repeating waveform, you can most easily find the constituent component parameters by aligning the Discrete Fourier Transform (DFT) frame on a whole number of cycles.  Suppose you chose four.  Then bins 0, 4, 8, 12, etc. contain the relevant values.

For a signal:

$$ x[n] = A \cos( \omega n + \phi ) $$

If

$$ \omega = \frac{2\pi}{N} k $$

where $k$ (cycles per frame) is an integer, then the $k$th bin will have a value of

$$ X[k] = \frac{A}{2} e^{i\phi} $$

So, you can read the tone's parameters right from the magnitude and angle of the DFT bin.  (FFT is an implementation of DFT)

If you aren't aligned on a whole number of bins, there is an interative process you can do with the formulas found in my blog articles.  This takes considerably more calculations though.


How do you construct common waveforms (square, sawtooth, triangle, noise) from sine wave partials using an IFFT? In the lovely image below we see the addition of sine waves making a square wave. I am coding a wave table synthesiser and thought it would be great to get the user to be able to fiddle with the values of the sine partials and then run an inverse FFT on them to produce the wave form.

I have some code for an FFT and you specify the number of buckets, real and imaginary arrays and the array length. The way it is used in the example is:

Imaginary array has all it's values set to 0
The real array has its 0th element set to 0
The real array has its arrayLength/2 element set to 0
Other values are set in the real array. Every time a value is set in the first half of the array one is set at the opposite end of the array, mirroring whatever is happening in the first half of the real array.
This spits out a single cycle of the waveform constructed using these sine partials

What are some common patterns for filling these partials to generate waveforms. How do you generate a square wave, a triangle wave, a sawtooth? Can you create noise or is that impossible if it is non periodic?

Extra: What are some tricks for producing 'good' waveforms? This is an opinion question, just a bonus if you have something to share about producing waveforms that sound good. (Good in your opinion :)

For those who would find it interesting there is a flash app here that lets you set the levels of the partials in order to construct waveforms, it's fun.

You can find the formulas for a rectangular, square, triangle, sawtooth and rectified sinusoidal waveform in Chapter 13 section The Fourier Series, of the Scientist and Engineer's Guide to Digital Signal Processing by Steven W. Smith.

If you IFFT a vector of random numbers, you get out another vector of random numbers that approximately follow a Gaussian distribution, due to the central limit theorem and that each time domain number was formed as a linear combination of multiple frequency domain numbers. But it will be periodic so it won't sound like pure noise unless the vector length is huge.

If you want to construct your own arbitrary waveform from line segments you can use:
$$f(x) = \left\{\begin{array}{ll}a + b\ x & \text{if}\ x_0 <= x < x_1\\ 0 & \text{otherwise}\end{array}\right.\\
\hat f(\omega) = \int^{+\infty}_{-\infty} e^{-i\omega\ x} f(x) = \int^{x_1}_{x_0} e^{-i\omega\ x}(a + b\ x) = -\frac{b\ \cos(\omega\ x_0)}{\omega^2} - \frac{(a + b\ x_0)\ \sin(\omega\ x_0)}{\omega} + \frac{b\ \cos(\omega\ x_1)}{\omega^2} + \frac{(a + b\ x_1)\ \sin(\omega\ x_1)}{\omega} + i\ \left(-\frac{(a + b\ x_0)\ \cos(\omega\ x_0)}{\omega} + \frac{b\ \sin(\omega\ x_0)}{\omega^2} + \frac{(a + b\ x_1)\ \cos(\omega\ x_1)}{\omega} - \frac{b\ \sin(\omega\ x_1)}{\omega^2}\right),$$

where $a = \frac{x_0\ y_1 - x_1\ y_0}{x_0 - x_1}$ and $b =\frac{y_0 - y_1}{x_0 - x_1}$ for a line segment starting at $(x_0, y_0)$ and ending at $(x_1, y_1)$ and  $\hat f(\omega)$ is the Fourier transform of $f(x)$. Discretely sampling $\hat f(\omega)$ at multiples of $2\pi$ will make $f(x)$ periodic by 1. If you have more than one segment you can sum their Fourier transforms.


FIR Filter 3dB Frequency I am new to signal processing and I want to calculate the 3 dB - cut off frequency of this FIR filter, which has the following transfer function:

For the evaluation i need to replace the z with  and set the squared magnitude equal to 1/2.

But the problem is when I want to solve this for omega zero is the result. My question is what am I doing wrong here, or is something missing ?

There are several problems with your approach:

The expression $\displaystyle\frac{1+z^{-1}}{2}$ is the transfer function $H(z)$, NOT its magnitude $|H(z)|$.
Consequently, the expression $\displaystyle \frac{1+e^{-j\Omega}}{2}$ is the complex frequency response $H(e^{j\Omega})$, not its magnitude.

So now you have to compute the squared magnitude of $H(e^{j\Omega})$ and set it equal to $\frac12$. This is a basic exercise in complex numbers, so I trust that you can take it from here.
For some magical reason, the right half of you last equation is actually the final correct equation, even though everything before it is wrong. If you think about it for a while, you should figure out that the solution is given by
$$\cos(\Omega_c)=0\tag{1}$$
which certainly does not imply that $\Omega_c=0$. What is it that it does imply?


Adaptive equalizer for IIR channel? I am quite new to the idea of equalization. I have a few queries regarding the same.

In my application I require to equalize a channel whose impulse response is an IIR response. Is it possible to design an adaptive equalizer based on LMS algorithm to equalize it? I tried some MATLAB simulations and it worked well (the estimation error reduced with each step and converged) for FIR channels but not for IIR responses. So is there any restrictions on the channel characteristics (like FIR/IIR response) for the adaptive filter to function properly?
Also, from my simulations for FIR channels I observed that, in absence of noise, the adaptive filter after converging had a response that is inverse of the channel response which is the same as that of a Zero-Forcing equalizer. So is it always true that, LMS adaptive filter will tend to a Zero-forcing equalizing filter in the absence of noise? If yes, can someone refer me to the proof of this?

I have another doubt. Say I am measuring the channel output in presence of an additive zero mean noise,

$$Y(z)=X(z)H_c(z)+N(z),$$

Where $H_c(z)$ is the channel response, I still would like the equalizer to estimate the inverse response of the channel. I thought this should work, as if I try to derive the Weiner-Hopf filter, in presence of noise, I end up getting

$$R_{dx}(k) - R_xW -E\left[n(l)x^{*}(n-l]\right],$$

Where $R_{dx}$ is the cross correlation between the desired signal $d(n)$ and the filter output $x(n)$, and $R_x$ is the auto-correlation matrix, and $W$ is the filter weights. Since the noise is independent of $x$ and has zero mean, $E\left[n(l)x^{*}(n-l)\right]$ becomes 0 and I again get the same equation. So I thought that the same algorithm could potentially work.

Also in the LMS algorithm, instead of using $E\left[e(n)x^{*}(n)\right]$, directly as $e(n)x^{*}(n)$, I tried to take an $L$ point average of the product, with the intention that, since the error $e(n)$ will consist components both due to ISI as well as the additive noise, averaging it will remove the additive noise and retain the ISI, and so the filter coefficients will converge to give the inverse response of the channel.

But MATLAB simulations do not agree with my argument above. I want to know the mistake in my argument and also I would like to know if there is any means by which I could still estimate the inverse response of the channel to a good approximation even in presence of a zero mean noise.

Note that the inverse of an FIR system is IIR, and the same is true for the inverse of an IIR system, unless it is an all-pole system, the inverse of which would be FIR. So in most cases the ideal equalizer should have an infinitely long impulse response in order to perfectly invert the channel. In practice almost all adaptive equalizers are FIR filters which can only partly compensate for the channel response. This is why I think that there should be no reason for your equalizer to perform worse for IIR channels than for FIR channels. You could post some more information about the channel responses and the convergence of your equalizer to make it easier to see what's going on.

As for the zero forcing criterion, such criteria for designing an equalizer only make sense for noisy channels, because otherwise there is no trade-off between equalizing the channel response and noise enhancement. Zero-forcing and MSE criteria result in the same equalizer in the absence of noise. If you have no noise, the only problem is the (hopefully linear) distortion of the channel and naturally this is the only thing the equalizer would compensate for. However, with an FIR equalizer exact zero-forcing is usually impossible as discussed above.


Why does the general formula for DFT of real sinusoids not apply here? In this course (Coursera: Audio Signal Processing for Music Applications) the professor derives a general equation for finding the DFT for real sinusoids:

  DFT of real sinusoids
  \begin{align}
x_3[n]&= A_0\cos\left(2\pi k_0 n/N\right)=\boxed{\frac{A_0}{2}e^{j2\pi k_0 n/N} +\frac{A_0}{2}e^{-j2\pi k_0 n/N}}\\
X_3[k]&=\sum_{n=-N/2}^{N/2-1} x_3[n]e^{-j2\pi k n/N}\\
&=\sum_{n=-N/2}^{N/2-1}\left(\frac{A_0}{2}e^{j2\pi k_0 n/N} +\frac{A_0}{2}e^{-j2\pi k_0 n/N}\right)e^{-j2\pi k n/N}\\
&=\sum_{n=-N/2}^{N/2-1}\frac{A_0}{2}e^{j2\pi k_0 n/N}e^{-j2\pi k n/N} +\sum_{n=-N/2}^{N/2-1}\frac{A_0}{2}e^{-j2\pi k_0 n/N}e^{-j2\pi k n/N}\\
&=\sum_{n=-N/2}^{N/2-1}\frac{A_0}{2}e^{-j2\pi\left(k-k_0\right) n/N}+\sum_{n=-N/2}^{N/2-1}\frac{A_0}{2}e^{-j2\pi\left(k+k_0\right) n/N}\\
&=N\frac{A_0}{2}\text{ for }k=k_0,-k_0; 0\text{ for the rest of }k
\end{align}

Keeping this in mind, given a discrete signal $x[n]$ that has 4 samples given as: $[1, -1, 1, -1]$

the equation from the slide can be used to represent our signal as:
\begin{align}
x[n] &= {A_0\over2}e^{j2 \pi k_0 n/N} + {A_0\over2}e^{-j2 \pi k_0 n/N}\\
&= {e^{j 4 \pi  n/4}\over2} +  {e^{-j 4 \pi  n/4}\over2}\\
&= {e^{j \pi  n}\over2} +  {e^{-j \pi  n}\over2}
\end{align}

And now again, calculating $X[2]$ with the equations from the same slide:

\begin{align}
X[2]&=\sum_{n=-N/2}^{N/2-1} \left({e^{j \pi  n}\over2} +  {e^{-j \pi  n}\over2} \right)  e^{-j2\pi kn/N}\\
&=\sum_{n=-N/2}^{N/2-1} \left({e^{j \pi  n}\over2} +  {e^{-j \pi  n}\over2} \right)  e^{-j4\pi n/4}\\
&=\sum_{n=-N/2}^{N/2-1} \left({e^{j \pi  n}\over2} +  {e^{-j \pi  n}\over2} \right)  e^{-j\pi n}\\
&=\sum_{n=-N/2}^{N/2-1} {e^{j (\pi n - \pi n)}\over2} +  {e^{-j 2 \pi  n}\over2}\\
&= \frac N2 + \frac N2\\
&= N
\end{align}

But this contradicts the last line on the slide that says $X[2]$ should be $N{A_0\over2}$.

Are the slides wrong?
I'm also having trouble figuring proving/understanding the following equality from the lecture slide:

  $$\sum_{n=-N/2}^{N/2-1} {A_0\over2} e^{-j2\pi(k-k_0)n/N} + \sum_{n=-N/2}^{N/2-1} {A_0\over2} e^{-j2\pi(k+k_0)n/N} = N\frac{A_0}{2}\text{ for } k=k_0, -k_0$$

PS: From the same course, some slide from a another lecture show these equations that are consistent with my derivation (see equation for $<x,s_2>$)

  DFT: scalar product:
  $$
\boxed{<x, s_k>=\sum_{n=0}^{N-1}x[n]s_k^{*}[n]=\sum_{n=0}^{N-1}x[n]e^{-j2 \pi k n/N}}
$$
  Example:

  $$x[n]=[1, -1, 1, -1]; N=4$$

  \begin{align}
<x, s_0>&=1*1 +(-1)*1 + 1*1 + (-1)*1 =0\\
<x, s_1>&=1*1 +(-1)*(-j) + 1*(-1) + (-1)*j =0\\
<x, s_2>&=1*1 +(-1)*(-1) + 1*1 + (-1)*(-1) =4\\
<x, s_3>&=1*1 +(-1)*j + 1*(-1) + (-1)*(-j) =0\\
\end{align}

It's not a contradiction. It's just that for an alternating signal the two contributions at $\pm k_0$ coincide and add up. Note that in your example with $N=4$, you have $k_0=\pm 2$, but these two values of $k_0$ are in fact identical due to the periodicity of $N=4$: $X[-2]=X[2]$.

The same would be true for $k_0=0$, i.e. a constant signal (with frequency zero and amplitude $A_0$). In that case you'd get $X[0]=A_0N$, just like in the case of a signal at Nyquist (maximum frequency, as in your example) where you get $X[N/2]=A_0N$. For all other values of $k_0$ you get two contributions with amplitudes $A_0N/2$.

Regarding your second question, this is an identity that is very important in DSP. First of all note that due to periodicity of the complex exponential function you have the equality

$$\sum_{n=-N/2}^{N/2-1}e^{-j2\pi (k-k_0)n/N}=\sum_{n=0}^{N-1}e^{-j2\pi (k-k_0)n/N}\tag{1}$$

For $k=k_0$ the complex exponential in $(1)$ becomes $1$, so the sum equals $N\cdot 1=N$. For $k\neq k_0$ you can use the formula for the geometric series to see that

$$\sum_{n=0}^{N-1}e^{-j2\pi (k-k_0)n/N}=\frac{1-e^{-j2\pi (k-k_0)}}{1-e^{-j2\pi (k-k_0)/N}}=0\tag{2}$$

because for $k\neq k_0$ we have $e^{-j2\pi (k-k_0)}=1$.


Combiner vs mixer vs coupler? All these three terms are often used in communication texts and apparently they are some how similar
But what is exact difference between each of them?

These are typical references to RF analog components where the three of these terms are used in the same context. They are often passive devices but can be implemented as active devices (requiring external power) as well. In the passive RF device world they are the following:

Combiner: This RF device combines two or more signals to a single output and is further described as a "Power Combiner". It is typically a passive reciprocal device: it can be used in the other direction as a "Power Splitter".  This can be implemented as a resistive combiner/splitter (for example put three 50 ohm resistors in "delta form" to create a resistive two way power combiner/splitter which will have a 50 ohm match when loaded at all ports, 6 dB loss between each port and no isolation. Often they are implemented with transformers at lower RF frequencies (10MHz to 2GHz), or with 1/4 wave transmission lines at higher frequencies (Wilkinson combiner/splitter) with internal isolation resistor(s) such that they have isolation between the input ports.

Mixer: This RF device multiplies the time domain signal presented at two ports and provides the product out. It typically will have an RF, IF and LO port corresponding to "Radio Frequency", "Intermediate Frequency" and "Local Oscillator Frequency". It is often used in frequency translation and phase detection applications as well as many other uses. This is often implemented as a passive device with Schottky diodes in a ring with balanced transformers as a "double-balanced mixer" but has many other forms. As a passive device it is reciprocal (any port can be an input or output) but there are often different frequency ranges for each port.

Coupler: This RF device has an input port, output port, and coupled port. The signal presented to the input port appears at the output port with minimal loss, while some of that signal is "coupled" off and appears at the coupled port.  The 2-way power divider (which is the 2-way power combiner run backwards) is a 3dB coupler and couples off 1/2 the signal power (although that is exactly what it is, we would typically not call it that). More common couplers are -6 dB, -10 dB, and -20 dB. The theoretical loss of a -10 dB coupler from input port to output port would be -0.4 dB, corresponding to the -10 dB signal that got coupled off. Further, the coupled port is isolated from the output port--- a signal presented there, such as a reflection of the coupled signal or otherwise, would be greatly attenuated (~ 30 dB) at the output port but would travel back to the input port. Therefore you can use a coupler in reverse to inject a signal into a system where the signal is fed into the output, the coupled input is fed into the coupled port and the combined signal will come from the input port with the relative coupling ratio. Similar to the power combiner, this can be implemented with transformers for frequencies from 10MHz to 1 GHz or adjacent transmission line structures at higher frequencies (strip-line and micro-strip forms).  There are also 4 port “bidirectional couplers” which have a forward and reversed coupled port. Basically the 3 port devices simply have the reverse coupled port internally terminated, but many of the bidirectional couplers are made as two 3-port couplers back to back for higher isolation performance.


Bits - symbols mapping in M-PAM beginner question I am trying to implement a M-PAM modulator in Matlab without using the built-in functions. I am confused when it comes to the symbols-bits mapping.

Say I have a sequence of bits: 1000110110 and M = 4 which means I encode k = log2(M) bits into each symbol. So, I get the vector s = [2 0 3 1 2]. But to my understanding this is not the proper mapping because the symbols arent symmetrical e.g. 1,-1,2,-2,3,-3. However, in functions like pammod I see examples that produce the symbols with randi that gives the first kind of mapping. Then, the mapping is passed as an argument in pammod. That is really odd for me since later on the amplitude per symbol is computed using the symmetrical mapping like (2m - M - 1)A (2m - M - 1 produces symmetrical symbols) where A is an energy constant and m = 1,2,3...,M.

So is it ok to use the first mapping when computing the Am amplitude but its less viable? Or I am misunderstanding the concept?

$2^k$-ary PAM maps $k$ bits (which can be expressed as integers in the range  $[0, 2^k-1]$ to $2^k$  levels that are equally-spaced and symmetric about $0$. ($0$ is NOT one of the levels).. The levels are odd integers, and so for $k=2$, the levels are $-3, -1, +1, +3$ which are equally spaced and symmetric about $0$.  As to which integer in $[0,2^k-1]$, there are several common conventions such as:
\begin{align}
00 &\to -3\\
01 &\to -1\\
11 &\to +1\\
10 &\to +3
\end{align}
or
\begin{align}
00 &\to -3\\
10 &\to -1\\
11 &\to +1\\
01 &\to +3
\end{align}
or
\begin{align}
10 &\to -3\\
11 &\to -1\\
01 &\to +1\\
00 &\to +3
\end{align}
etc.
In all cases, there is a one-bit change in going from one level to the next. Note also that your formula $2m-M-1, m = 1,2,\ldots, M$ corresponds to the first list if you assume that $m-1$ is expressed as a $\log_2M$-bit integer and expressed in Gray code so that adjacent symbols differ in only one of the $\log_2M$ bits.


Question about z transform After studying z transform from different books and literature on internet I want to ask few which makes me confuse.

a) From the Discrete Time Fourier Transform we have drive equation for z transform. $$ X(z)= \sum _ {n=-\infty}^{+\infty} x[n]z^{-n}$$ where $z$ is represented in polar form $z=re^{j\omega}$
I want to know that why we represent $z$ in polar coordinates? as in some books it is written than $z$ is complex $z=\sigma + j \omega$

b) ROC for $z$ transform is same as Laplace Transform? In laplace transform we check that direction of $t$ (i.e. if we have $u(t)$ than the $Re[s]  > a$)?

The variable $z$ is complex, and it can be represented either in polar form $z=re^{j\phi}$ or in terms of its real and imaginary parts. The polar form is preferred because the region of convergence (ROC) is determined by $r$, the magnitude of $z$. This becomes clear if you write the term inside the sum of the definition of the $\mathcal{Z}$-transform as

$$x[n]z^{-n}=x[n]r^{-n}e^{-jn\phi}\tag{1}$$

From $(1)$ it is clear that the value of $r$ determines the convergence of the series. This is why the ROC of the $\mathcal{Z}$-transform is specified in terms of $r=|z|$. For the bilateral $\mathcal{Z}$-transform the ROC is a ring ($a<|z|<b$, with $b>a\ge 0$), whereas for the unilateral $\mathcal{Z}$-transform, the ROC is outside a circle centered in the origin of the complex plane ($|z|>a$, $a>0$).

Note that for the Laplace transform the complex variable $s$ occurs in the exponent, so it's the real part of $s$ that determines the ROC of the Laplace transform. These are the correspondences between the Laplace transform and the $\mathcal{Z}$-transform in terms of ROCs:

$$\begin{align}&s\text{-plane }&&z\text{-plane}\\\\
&\text{right half-plane} && \text{region outside a circle (centered at the origin)}\\
&\text{vertical strip} && \text{ring (centered at the origin)}\\
&\text{left half-plane} && \text{region inside a circle (centered at the origin)}
\end{align}$$

The $j\omega$-axis in the $s$-plane corresponds to the unit circle in the $z$-plane, i.e. if the unit circle is part of the ROC, then the DTFT exists, and if the corresponding sequence is the impulse response of an LTI system, then the system is stable.


BER Analysis of Dual hop Decode and Forward Relay I want to perform fair BER analysis between Two-hop Decode and Forward relay communication (S-->R-->D) and Direct Transmission from S-->D. The links have Rayleigh fading.
Now for fair analysis 1) Total Transmit Power is kept same 2) Rayleigh Channels with  distance are considered for each link.

In the literature,BER vs (Average SNR) is plotted.

How do I calculate Average SNR for dual hop transmission?

and Can I plot BER vs Transmit power?

I am working with a similar kind of network. I don't know if it is correct, but I think the most proper way is to consider the SNR from the first transmitter, not considering what the relay do to the signal. This because you want to analyze the performance of the end-to-end communication. SO you send the signal with some energy, and you want to know the performance of your signal once you have sent it with that energy.

That's the way I am doing, and comparing it to the direct transmission, I am getting a very similar performance: a difference of less than 1dB for the same BER.

I am not amplifying the signal in the relay, so I am just doing the decode and forward (DF). I think it is reasonable to get this "worst" performance in the DF relay compared to the direct transmission because you have two channels, therefore you have two components of noise. You perform the DF in the relay, but it can decode wrong.

Let me know if this answer was useful, and if you have any more questions. I am not an expert in this area, but I am working on it.


Jacobian of affine warp in Lucas-Kanade image alignment I hate this kind of questions, but I'm really stuck. I'm trying to implement Lucas-Kanade algorithm as described in this paper (see pages 4 and 5). Unlike most explanations I've seen they don't assume optical flow, but instead use affine transformation for most examples and in general refer to warp function just as $W(x, p)$, where $x$ is the point and $p = (p_1, p_2, p_3, p_4, p_5, p_6)^T$ is parameter vector of affine transformation.

At the beginning of page 4 authors outline their version of Lucas-Kanade algorithm. I'm stuck at steps (4) and (5), namely, evaluating the Jacobian $\frac{\partial W}{\partial p}$ and calculating the steepest descent images $\nabla I\frac{\partial W}{\partial p}$.

(As far as I can understand, $\nabla I$ is just a pair of 2 matrices $(\frac{\partial I}{\partial x}, \frac{\partial I}{\partial y})$ in this context. Please, correct me if I'm wrong).

Authors carefully describe affine warp and even provide formula for its Jacobian (Equation 8):

$$\frac{\partial W}{\partial p} = \pmatrix{x & 0 & y & 0 & 1 & 0 \\ 0 & x & 0 & y & 0 & 1}$$

However, Jacobian of the warp is (not surprisingly) defined only for a single pixel and not the entire image. But in step (5) we calculate multiplication of image gradient and this Jacobian - $\nabla I \frac{\partial W}{\partial p}$, and as far as I can see from the context (see, for example, Figure 2 on page 5), it is done for the whole image and not per pixel.

My question is, how should I interpret this multiplication and what are the real sizes/formats of $\nabla I$ and $\frac{\partial W}{\partial p}$?

I understand that this question may require reading a lot from that paper, so I will be glad to explain any point that is not clear enough. Also feel free to refine question title or contents if you have an idea of a better wording.

You're right that all the quantities are computed for a single pixel, so in the product

$$\nabla I\cdot \frac{\partial W}{\partial p}\cdot\Delta p\tag{1}$$

the sizes of the vectors and matrices are

$$\nabla I\quad\;\;\; 1\times 2\\
\frac{\partial W}{\partial p}\quad 2\times n\\
\Delta p\quad\;\;\;\; n\times 1$$

where $n$ is the number of parameters. So the expression in (1) is a scalar. The total error measure given by Equation (6) in the paper is the sum of this scalar expression over all pixels.


Compressing A Range Of Decibel Values to One "Representative" Value Don't have any formal dsp experience - sophomore in college.

Working on constructing a 25 band bark spectrogram for psychoacoustic analysis of audio FFT data. To do this, I'm grouping all FFT frequency bins into their relevant bark band and then averaging (summing them and dividing by the total number) the db values in each band to get the representative value for each band.

Is this the correct way to represent a bunch of db values as a single value?

The answer provided by Jason R in the comments --

  "dB values" are really just a way of representing values on a logarithmic scale. Arithmetic-averaging values on a logarithmic scale is equivalent to geometrically averaging the original values (in linear scale). This probably isn't what you want. You probably would like a measure of total power in each band, which you would get by just summing the power spectrum bins across each band. You could divide by the number of bins, but if the bands aren't equal-width, that will add some bias to the result.

does a good job of identifying the problem but I'm not sure it is a correct solution (or maybe I'm just reading it wrong). If a summation is to occur, it should occur on the raw magnitudes of the bins (before conversion to dB).

The issue is this

3dB + 3dB != 6dB

to do a summation on the decibel scale is more complex than simply adding the values together. It can be done but if so you should do some reading on the dB scale and the nature of adding relative power. Doing the addition when still in a linear domain is easiest.


What is the difference between a'mixer' and a ' multiplier' used in modulation process? I am not sure if it is valid to conclude that a mixer produces a difference/sum of frequencies at output ,while multiplier produces multiples of frequency at its output. I am keen to know how far I have got it right. Is there something else that I should know? please post with simple e.g.

For most purposes a mixer is a multiplier.

The issue regarding sum / difference frequencies is that when you multiply two signals of different frequencies, $\omega_1$ and $\omega_2$ say, then you get one component at the sum of the two frequencies and one component at the difference of the two frequencies:
$$
\cos(\omega_1 t) \cdot \cos(\omega_2 t) = \frac{1}{2} \left\{ \cos([\omega_1+\omega_2]t) + \cos([\omega_1 - \omega_2]t) \right \}$$

Response to comment

A frequency multiplier is a completely different kettle of fish. These are often highly nonlinear systems that take a signal $\cos(\omega_1 t)$ and produce $$\sum_{k=1}^{N} \cos(k \omega_1 t + \phi_k)$$

For example, a frequency tripler can be had by driving a transformer into saturation (so that it outputs something like a square wave) and low-pass filtering the result.


Does a raised-cosine filter have infinite support? I recently asked a question here regarding limiting the bandwidth of a digital baseband signal and was introduced to pulse-shaping using a "raised-cosine filter". At the time, I mistakenly assumed this to be a sinc function combined with the Hann window. This is because the Hann window is also known as the raised cosine window, and that seemed logical to me. However, after reading on Wikipedia, I see that the "raised-cosine filter" is actually another beast entirely. This function appears to be a sinc function windowed by another sinc-like function, with two parameters, β and T. The extent of this function, like sinc, appears to be infinite, for all values of β and T. This thoroughly confuses me, as I thought the whole point of windowing the sinc function was to give it finite support. What am I missing?

An ideal raised cosine impulse is band-limited with cut-off frequency $f_c=(1+\beta)/(2T)$, where $\beta$ ($0\le\beta\le 1$) is the roll-off factor. Since it is band-limited, it must have infinite support in the time domain. The important difference with a sinc impulse is that its envelope decays much faster (for $\beta>0$, of course).

Note that even though an ideal raised cosine filter has infinite support, in a practical implementation you would use a filter with finite support (implemented as an FIR filter) that approximates the ideal raised cosine impulse in some interval around its maximum.


pi/4 QPSK modulation/demodulation in GNU Radio Is there any quick mechanism to create pi/4 QPSK or pi/2 BPSK modulation scheme in GNU Radio using Constellation Modulator object or Constellation Object? Or I need to write it by myself? Thanks in advance

use chunks to symbols to convert your data to PSK symbols, and a vector source to alternatively generate $0, \Delta \varphi$. Multiply the two streams, then do pulse shaping.


How to use the iircomb filter in matlab I have a sound file that needs filtering which looks like this:

Hence I think that I have noise spaced equally at 832Hz. i tried to use the iircomb filter taken from mathworks website and i managed to create this:

clear
fs = 44100;
fo = fs/53;
q = 35;
bw = (fo/(fs/2))/q;
[b,a] = iircomb(fs/fo,bw,'notch'); % Note type flag 'notch'
fvtool(b,a);
[f,fs] = audioread('sample3b_va18535.wav');
fnew=filter(b,a,f)

N = size(f,1);
df = fs / N;
w = (-(N/2):(N/2)-1)*df;
y = fft(fnew(:,1), N) / N; % For normalizing, but not needed for our analysis
y2 = fftshift(y);
figure(2);
plot(w,abs(y2));
p = audioplayer(fnew,fs);
p.play;

Now, the filter works pretty good but there is still a little bit of noise after filtering which can be heard. Should I just live with it, or perhaps my filter could be improved? Could someone with more experience work on it for a minute or two? The sound file is only 3 seconds long.

Here is the audio file if anyone needs it: SOUND FILE

Thank you for help.

Precision delay is, in my opinion, best understood in terms of the Nyquist-Shannon sampling and reconstruction theorem.

If a continuous-time (a.k.a. "analog", but that is somewhat imprecise) signal $x(t)$ is bandlimited to below the Nyquist frequency $\frac{1}{2T}$ (where $T$ is the sampling period and $f_\mathrm{s} \triangleq \frac1T$ is the sample rate), then the original continuous-time signal $x(t)$ can be reconstructed from the discrete samples $x[n]$ as:

$$\begin{align}
x(t) &= \sum_{n=-\infty}^{\infty} x(nT) \, \operatorname{sinc}\left(\tfrac{t - nT}{T}\right) \\
&= \sum_{n=-\infty}^{\infty} x[n] \, \operatorname{sinc}\left(\tfrac{t - nT}{T}\right) \\
\end{align}$$

where

$$ \operatorname{sinc}(u) \triangleq \begin{cases}
\frac{\sin(\pi u)}{\pi u},  & \text{if } u \ne 0 \\
1, & \text{if } u = 0 \\
\end{cases} $$

All terms of the summation are bandlimited to a maximum frequency of $\frac{1}{2T}$, so the summation is bandlimited to the same bandlimit as the original $x(t)$.  The samples are

$$ x[n] \triangleq x(t) \Bigg|_{t = nT} \triangleq x(nT)$$

Because $x(t)$ can be evaluated at any $t$, not only multiples of $T$, this is the basis in which we interpolate to get a precision delay.

Suppose we wanna delay by an amount of time $\tau$, then the output of the delay is

$$\begin{align}
y(t) &= x(t-\tau) \\
&= \sum_{n=-\infty}^{\infty} x[n] \, \operatorname{sinc}\left(\tfrac{t-\tau - nT}{T}\right) \\
\end{align}$$

So that is explicitly where you get the delayed signal of $y(t) = x(t-\tau)$ from the samples $x[n]$.  The first problem we have is that we cannot add up an infinite number of terms.  So our first practical approximation is truncating terms.  But that is the same as applying a rectangular window (usually considered to be the worst one), so instead, we'll apply a better window.

$$\begin{align}
y(t) &= x(t-\tau) \\
&= \sum_{n=-\infty}^{\infty} x[n] \, \operatorname{sinc}\left(\tfrac{t-\tau - nT}{T}\right) \\
&= \sum_{n=-\infty}^{\infty} x[n-n_\tau] \, \operatorname{sinc}\left(\tfrac{t-\tau - (n-n_\tau)T}{T}\right) \\
&= \sum_{n=-\infty}^{\infty} x[n-n_\tau] \, \operatorname{sinc}\left(\tfrac{t-(\tau-n_\tau T)}{T}-n\right) \\
\end{align}$$

... for any integer $n_\tau$.  Let's choose it as the integer part of the dimensionless $\frac{\tau}{T}$, which is the delay in sample units.

$$ n_\tau \triangleq \left\lfloor \frac{\tau}{T} \right\rfloor $$

The fractional part of the delay (in sample units) is

$$ \frac{\tau}{T} - \left\lfloor \frac{\tau}{T} \right\rfloor $$

(i gotta return to this a little bit later.)


Is it possible to design noncausal filters in Matlab? I'm trying to design a lowpass IIR filter with 5Hz cutoff frequency (sample frequency is 200Hz). i know how to create it with fdatool, or with functions such as 'butter' or 'ellip'. However, i think that all the filters created are causal. Is it possible to create also noncausal filters? Thanks in advance

The short answer is: yes, of course it's possible to design non-causal filters. Consider a simple first order IIR filter:

$$y[n]=ay[n-1]+bx[n]\tag{1}$$

If implemented in the way as indicated in Eq. (1), this filter is causal. It is also stable if $|a|<1$. Now rewrite (1) as

$$y[n-1]=\frac{1}{a}y[n]-\frac{b}{a}x[n]\tag{2}$$

Implemented in this way, you compute $y[n-1]$ from the future samples $y[n]$ and $x[n]$, so the system is anti-causal. Note that it is stable if $|a|>1$. Similar considerations apply to higher-order filters.

One way of realizing non-causal IIR filters is implemented in Matlab's filtfilt function. How this works is described in this answer.

Non-causal FIR filters can always be implemented by appropriate buffering of the input signal and delaying the output signal.


Finding filaments in high dimensional space I have unit vectors in a high-dimensional space (~300) that are non-uniformly distributed over this hypersphere. There is physical meaning when the vectors are clustered together, but I'm hypothesizing that there is additional meaning when the vectors form a "path" along this hypersphere.

For the purposes of this question, let's assume that the deviation between the geodesic across the hypersphere and a straight line are are approximately equal (that is, I'm assuming the distances are small relative to the radius).

The question becomes, "From a large 300-dimensional point cloud in Euclidian space, how do I detect groups of points that are aligned in any direction?"

Following my comment, I now have time to post it as an answer.
If you expect that a large enough portion of the vectors belong to that line, you can use RANSAC.
RANSAC works by selecting a small and random subset of the data and fit it to the desired model (say a line). The fitting can be achieved in any way you like, e.g. least-squares. Then, the data is divided to inliers (the points in space that "agree" with the model) and outliers (those who don't "agree"). The process is repeated until you have enough inliers.
This would only work if the probability of selecting the points on the line is high enough.

As far as multiple lines, I don't think RANSAC supports this inherently. You can do it iteratively by removing all inliers from the set once a line has been found.

If you are working in MATLAB\Octave, my favourite RANSAC implementation can be downloaded from here.


What is the difference between wide sense and strict sense stationary processes? What is the difference between wide sense and strict sense stationary processes (SP) ?

According to the definition (by Heinrich Meyr, Marc Moeneclaey, Stefan A. Fechtel in ￼"Synchronization, Channel Estimation, and Signal Processing") :

strict-sense SP = not time dependent
wide-sense SP = not dependent on variable t (time)

Could someone please illuminate what this really means ?

Perhaps pointing out a good book in which the difference is explained ?

If an SP is not dependent on variable t (time) then how can it depend on time ( a process which is not strict sense stationary but wide sense stationary)?

Could you give an example of such a process ?

Why and when is making this difference useful/neccessary?

A process is stationary if:

its mean is a constant value: $\mu_x(t)=\mu x$
its MSV(mean square value) is a constant value.
its variance is a constant value. $\sigma^2_x(t)=\sigma^2_x$
its autocorrelation depends on the time distance between 2 samples : $$R_x(t_1,t_2)=R_x(\tau)$$
its auto-covariance depends on the time distance between 2 samples at time $t_1$ and $t_2$ : $$K_x(t_1,t_2)=K_x(\tau)$$
So, with respect to information above, we cannot just say that a process is stationary if it does not depend on time, because for mean,MSV and variance a stationary process does not depend on time since its values are constant, but for autocorrelation and auto-covariance it does depend on time actually, but the function itself does not depend on time since it depends only on the time distance between 2 samples.

The process is SSS (strict-sense stationary) if its $N^{\rm th}$ order probability density function is stationary for any given $N$.

The process is WSS (wide/weak-sense-stationary) if its mean value is constant $\mu_x(t)=\mu_x$, and its autocorrelation function depends on time difference between 2 samples, $$R_x(t_1,t_2)=R_x(\tau).$$

Important:
Remember that all SSS processes are also WSS, but NOT all WSS processes are SSS.There is only one exception: "Jointly Gaussian Processes". Jointly Gaussian processes have any order of probability density functions(property of SSS). So, even if a jointly gaussian is WSS it is also SSS.And also remember that white gaussian noise process is a jointly gaussian process with 0 mean. So it is WSS and also SSS. Since its mean value is 0, the we can say that its autocorrelation function and autocovariance function are equal and they do depend on time difference between 2 samples  at time $t_1$ and $t_2$. $R(\tau)=K(\tau)$ where $\tau=t_1-t_2$.

In conclusion,

if you see only the mean value is constant (not depend on time), and the autocorrelation function depends on only time difference $\tau= t_1-t_1$. $R(\tau)$. --> WSS
if you see that the $N$-th order pdf of process is stationary for any $N$ --> SSS
if the process is jointly gaussian --> WSS and SSS
if the process is white gaussian noise process --> WSS and SSS with mean=0 and $$R(\tau)=K(\tau).$$


Binary robust independent elementary features ("BRIEF") understanding How are the descriptors in the BRIEF algorithm matched to each other?

How is the template image found in the another image? How do we compare the descriptors?
I have read that article but didn't understand how they have done it.

They wrote that they compared descriptors by using Hamming distance, but how is invariance to changes in illumination, size, etc. achieved?

First of all, there's no such thing as 'template' in this paper - the word 'template(s)' has a different meaning in Computer Vision.

The method used in this paper is relatively straight-forward. Let me break it down for you. There are three important things that you need to do when doing tasks such as object recognition, image matching, image stitching, and similar others, using Local Features

The first important thing is Detection; in this step you want to detect points-of-interest or keypoints, and what that means is that you want to choose local points (basically small patches) that you think are interesting in the image, there are many ways to do that; this paper doesn't contribute in this area. However, it seems that they use SURF feature detector and CenSurE keypoints (you can look them up if you want to know more about them, I won't talk about this step except that they use features such as gradients and so, which means that if you have a photo, it's unlikely that a point in the middle of the sky would be chosen as an interest point, that's because the pixels around that point are all same intensities, things that are 'busy' tend to be chosen as interest points (e.g. table/building edges/corners)).

After detection is done, Feature Description follows. You know the interesting points in the image and now you want to describe them (basically you want to describe the points/patch around the interesting points). SIFT is one popular feature descriptor. In this paper, they come up with a new one, called BRIEF. BRIEF is based on comparisons, so let's say we have a patch (50 pixels by 50 pixels), we choose two points and compare the intensities of the two points, if the first point is larger than the second point, we assign the value '1', else '0', we do that for a number of pairs and we end up with a string of boolean values. Now the big question is 'how do you pick the pair of points?', in the paper they explained 5 ways, relatively similar, I will describe the first one. What you do is you uniformly (equal probability) choose a point from -S/2 to S/2, in our example we said that the patch size is 50, so we choose a point between -25 to 25. Assuming that the 0,0 coordinate is located at the center of the patch. So here's an example;

We want to select the first pair, each point consists of (X,Y) coordinates so we uniformly select the X-coordinate of the first point, then the Y-coordinate of the first point, let's assume it's (10,-1), now the second point; we uniformly select the X-coordinate for the second point, and the Y-coordinate for the second point, let's assume it's (-2,20), now we get the intensity values for each point and see which one has the larger intensity values - if first is larger we assign the first boolean value to be '1', if not, we assign '0'. We do that for many pairs, and we end up with a vector of boolean values.

*Very important point: I assume for this to work, you will have to specify a seed value before the random generator. What that means is that you want to choose the same values for each patch - this is super important because when you compare/match patches, the whole system will break down if they're not chosen the exact same way. - hopefully this will make sense when you read the matching step.

So, we do that for each interest point that was detected by the detector. For each interest point, we will get a vector of boolean values.

Now to match two images (third step; matching), we do the exact same thing for the other image, we detect, then describe using the BRIEF. For example let's say we have 10 interest points in each image (this can always work if we get the 10 most interesting points in each image), we use BRIEF to describe each patch using for example 50 pairs, so each image will be described by 10 vectors of 50 boolean values.

To compare the two images, we find the nearest neighbor of each vector from the first image to each other vector from the second image. We use Hamming distance which is pretty fast, example of hamming distance

hammingDistance((0, 1, 1), (0, 0, 0)) = 2

hammingDistance((0, 1, 1), (0, 1, 1)) = 0

hammingDistance((0, 1, 1), (1, 0, 0)) = 3

hammingDistance((0, 1, 1), (1, 1, 1)) = 1

Basically how many wrong correspondences

Hope that helps


Determining output of a LTI system Consider following LTI system $$y[n] - 2y[n-1] = x[n]$$ where $x[n]$ is the input to the system and $y[n]$ is the output. Let $x[n] = \cos[n\frac{\pi}{3}] + 2\cos[n\frac{\pi}{2} + \frac{\pi}{4}]$, determine the Fourier series coefficients for $y[n]$.

My try: First of all we should determine impulse response $h[n]$. So we have $h[n] - 2h[n-1] = \delta[n]$ and assuming initial rest condition leads to $h[n] = 2^nu[n]$. Then the frequency response could be computed as follows:
$$H(e^{j\omega}) = \sum_{n = -\infty}^{+\infty}h[n]e^{-jwn} = \sum_{n = 0}^{+\infty}2^ne^{-jwn} = \sum_{n = 0}^{+\infty}(2e^{-jw})^n$$
It's a geometric series and divergent because $|2e^{-j\omega}| = 2 >1$. So we should conclude that response to $x[n]$ doesn't exist? Also is initial rest a valid assumption here?

Edit: Let $a_k$ be the Fourier series coefficients for the input signal $x[n]$. Then the Fourier series coefficients for output is $b_k = a_kH(jk\omega_0)$ where $\omega_0$ is the Fundamental frequency of the input. The problem is that for each $k$ the geometric series diverges so $b_k$ doesn't exist.

You are right that the system $h[n] = 2^n.u[n]$ is unstable in absolute summability sense. And absolute summability of impulse response of an LTI system is the required condition for BIBO stability. Meaning Bounded Input bounded output stability. Here the input is bounded and it is only corresponding to 4 frequencies in digital frequency domain $\omega = -\pi /3, \pi /3, -\pi/2, \pi/2$. The magnitude of input cannot be greater than 3 at any value of $n$. But, the output will not be bounded because the system is not BIBO stable.

Even if you start the system at 0, meaning $y[n] = 0, \forall n < 0$, $y[n]$ will blow up as $n$ grows.
$$y[n] = \sum^{\infty}_{m=0}2^m.x[n-m], \forall n \ge 0$$
Putting $x[n] = cos[\frac{\pi n}{3}]$ as input:
$$y[n] = \sum^{\infty}_{m=0}2^m. \frac{e^{j\pi (n-m)/3} + e^{-j\pi (n-m)/3}}{2}$$
$$y[n] = \frac{e^{j\pi n/3}}{2}.\sum^{\infty}_{m=0}(2.e^{-j\pi/3})^m + \frac{e^{-j\pi n/3}}{2}.\sum^{\infty}_{m=0}(2.e^{j\pi/3})^m$$
Since both the sums are going to blow up as $n$ grows, we can say that $y[n]$ is not bounded.

Another way of saying this is that $H(e^{j\omega})$ does not exist for this system ,because DTFT sum does not converge. And so we cannot get $|H(e^{j\omega})|$ and $\angle{H(e^{j\omega})}$ for any $\omega \in [-\pi, \pi]$.

If we considered the system to be anti-causal, then the unit circle will fall inside ROC because ROC will become $|z|<2$. And the anti-causal but stable system's impulse response will be given by : $h[n] = 2^n.u[-n]$.

But this will be wrong thing to do, because the difference equation represents a Causal System, since the current output sample depends only on current input sample and a past output sample.


Gain function calculation (frequency response) Define moving average process $y_t := 0.5 x_t + 0.5 x_{t-1}$ where $x_t := e^{i2 \pi t}$. Its frequency response is then:

$$H(f) = 0.5 + 0.5 e^{-i2\pi f}$$

Recall that the frequency response in polar notation is:

$$H(f) = G(f)e^{i \theta}$$

where $G(f)$ is the gain function. We then have that:

$$ H(f) = 0.5 \big( e^{-i \pi f} + e^{i \pi f} \big)e^{-i \pi f}$$.

QUESTION: Can someone please explain why the last step is true? An explanation with some detail is much appreciated (I am new to signal analysis).

One way to see the equivalence is to first factor out the $0.5$ in your initial transfer function, i.e.

$$
H(f) = 0.5(1 + e^{-i2\pi f}).
$$

Then you can rewrite the $1$ and $e^{-i2\pi f}$ as follows:

$$
H(f) = 0.5(e^{i\pi f}e^{-i\pi f} + e^{-i\pi f}e^{-i\pi f}) = 0.5(e^{i\pi f} + e^{-i\pi f})e^{-i\pi f}.
$$

Your result then follows as a trivial transform.


How to transform an aerial image This is the first time I'm diving a little deep into images and computer vision. I have an aerial image taken from a drone (I know the pitch, roll, yaw of the drone and the angle of the camera relative to north). How do I transform the image so that it looks as if it is taken from the top? I want to map latitude, longitude on the aerial image to pixels (I know the lat, lon of all four corners of the image) and my understanding is that I need to transform the image before I can interpolate and map a new lat,lng pair to pixels.

Suppose you're imaging a square with sides of length $d$, and you're a height $h \gg d$ above the square. You're tilted in one dimension by $\theta$. Then the center of the square appears displaced by $-h \sin \theta$ and its apparent width along that dimension is $d \cos \theta$, not $d$. So, what object would be imaged onto a single pixel of size $d_p$? It would need to be a rectangle, displaced from the image center by $h \sin \theta$. Its width along the tilted dimension would have to be $d / cos \theta$, but along the untilted dimension it should still be $d$.

So after correction, the image consists of rectangular pixels of width $d / \cos \theta$. There's no transform on the raw image data that will implement this directly. However, it's easy to write a function that will return the value of the image at a given lat/long.


Change the melody of human speech using FFT and polynomial interpolation I'm trying to do the following:

Extract the melody of me asking a question (word "Hey?" recorded to
wav) so I get a melody pattern that I can apply to any other
recorded/synthesized speech (basically how F0 changes in time).
Use polynomial interpolation (Lagrange?) so I get a function that describes the melody (approximately of course).
Apply the function to another recorded voice sample. (eg. word "Hey." so it's transformed to a question "Hey?", or transform the end of a sentence to sound like a question [eg. "Is it ok." => "Is it ok?"]). Voila, that's it.

What I have done? Where am I?
Firstly, I have dived into the math that stands behind the fft and signal processing (basics). I want to do it programatically so I decided to use python.
I performed the fft on the entire "Hey?" voice sample and got data in frequency domain (please don't mind y-axis units, I haven't normalized them)

So far so good. Then I decided to divide my signal into chunks so I get more clear frequency information - peaks and so on - this is a blind shot, me trying to grasp the idea of manipulating the frequency and analyzing the audio data. It gets me nowhere however, not in a direction I want, at least.

Now, if I took those peaks, got an interpolated function from them, and applied the function on another voice sample (a part of a voice sample, that is also ffted of course) and performed inversed fft I wouldn't get what I wanted, right?
I would only change the magnitude so it wouldn't affect the melody itself (I think so).
Then I used spec and pyin methods from librosa to extract the real F0-in-time - the melody of asking question "Hey?". And as we would expect, we can clearly see an increase in frequency value:

And a non-question statement looks like this - let's say it's moreless constant.

The same applies to a longer speech sample:

Now, I assume that I have blocks to build my algorithm/process but I still don't know how to assemble them beacause there are some blanks in my understanding of what's going on under the hood.
I consider that I need to find a way to map the F0-in-time curve from the spectrogram to the "pure" FFT data, get an interpolated function from it and then apply the function on another voice sample (I mean curve fitting)
**Is there any elegant (inelegant would be ok too) way to do this? I need to be pointed in a right direction because I can feel I'm close but I'm basically stuck.
I need to increase/manipulate frequencies at certain point in time, I know it does mean to move data from one bin to another but I need to know how to work with the STFT data.
**
The code that works behind the above charts is taken just from the librosa docs and other stackoverflow questions, it's just a draft/POC so please don't comment on style, if you could :)
fft in chunks:
import numpy as np
import matplotlib.pyplot as plt
from scipy.io import wavfile
import os

file = os.path.join("dir", "hej_n_nat.wav")
fs, signal = wavfile.read(file)

CHUNK = 1024

afft = np.abs(np.fft.fft(signal[0:CHUNK]))
freqs = np.linspace(0, fs, CHUNK)[0:int(fs / 2)]
spectrogram_chunk = freqs / np.amax(freqs * 1.0)

# Plot spectral analysis
plt.plot(freqs[0:250], afft[0:250])
plt.show()

spectrogram:
import librosa.display
import numpy as np
import matplotlib.pyplot as plt
import os

file = os.path.join("/path/to/dir", "hej_n_nat.wav")
y, sr = librosa.load(file, sr=44100)
f0, voiced_flag, voiced_probs = librosa.pyin(y, fmin=librosa.note_to_hz('C2'), fmax=librosa.note_to_hz('C7'))

times = librosa.times_like(f0)

D = librosa.amplitude_to_db(np.abs(librosa.stft(y)), ref=np.max)

fig, ax = plt.subplots()

img = librosa.display.specshow(D, x_axis='time', y_axis='log', ax=ax)

ax.set(title='pYIN fundamental frequency estimation')

fig.colorbar(img, ax=ax, format="%+2.f dB")

ax.plot(times, f0, label='f0', color='cyan', linewidth=2)

ax.legend(loc='upper right')
plt.show()

Hints, questions and comments much appreciated.

The low quefrency MFCC coefficients might be a place to start.
I read a few articles by Stylianou from a few years back where he did funky stuff to peoples voice (gender and age modification).