Efficient Marker Detection for Images I need to detect a known number of simple markers (chessboard patterns) in a color image efficiently.

I implemented normalized cross correlation (NCC) as a pixel shader via WebGL since I need this to work in a client application in a browser. Detection works fine and is sufficiently robust, but I need better performance. It currently takes a number of seconds to calculate the NCC image, and I'm already using 'obvious' optimizations such as cropping the image (I have a rough indication where the markers should be) and reducing the resolution. I also optimized the code itself (e.g. using loop unrolling, but that is another matter and probably belongs on SO).

The image has a resolution of 900x520 pixels and the markers are 22 square pixels. The code runs on tablet computers with Intel Atom processors.

My question is: Are there algorithms that are better at utilizing what I know about the image and the marker geometry, ideally ones that are suited for almost branchless pixel shaders? For instance, I was considering a band pass filter (a fourier transformation that only leaves the frequency of the marker, so to speak) or FNCC (JP Lewis's paper), or something that utilizes the fact that the markers are black and white.

Unfortunately, it seems these approaches are comparatively expensive since the FFT of the image isn't cheap and my markers are small compared to the image resolution. Testing a ton of algorithms blindly isn't an option, unless there's a library I'm not aware of that helps me implement this in JavaScript or as glsl pixel shaders. Even then, an indication of a more efficient algorithm would be helpful.

EDIT:
Added an image of the marker. The background is normally transparent, of course, i.e. the NCC will ignore those pixels

This is a saddle point. I would suggest using derivative operators to spot such discontinuities. In fact, Haralick presents a broad overview of these methods in his Topographic Primal Sketch.

The use of approximate Facet model would give you a speed boost. I used them for edge detection here.

There is also one implementation available here.

Even though this question and the topic is rather old, I would like to update that, the method presented in the following work contains no branching and is quite suitable for the shader:

  Birdal, Tolga, Ievgeniia Dobryden, and Slobodan Ilic. X-Tag: A
  Fiducial Tag for Flexible and Accurate Bundle Adjustment. 3D Vision (3DV), 2016 Fourth International Conference on. IEEE, 2016.


If a Kalman filter can only receive information on $(x, y)$ position, is there any reason to have acceleration as part of the model? I'm new to working with Kalman filters. I have a Kalman filter that's working well for modelling the 2D position of an object. The model at the moment holds information on position ($x$ and $y$) and velocity ($x$ and $y$). I only have position information available for measurement. I was wondering if adding acceleration to the model would improves things in any way, where the only measurement information available is still just the $x$ and $y$ position. Would adding acceleration to the model make it so sudden jerks in movement are handled better or would there be no difference at all?

Whether to include a particular parameter in the model depends on how much you know about the parameter. If you're already including 2D velocity in the state, then the acceleration will be coming in through the random noise assumptions you are making.

Whether tracking a person will improve if you include the 2D acceleration in the state depends on what assumptions you can make about how the person's acceleration changes.

For some models, we might know that the jerk$^1$ is completely random but not the acceleration. In that case, it would make sense to include the acceleration and feed the randomness in as the jerk (change to the acceleration).

${\tiny^1\sf \mbox{The derivative of acceleration, I'm not casting aspersions. }}$


How is ESPRIT pronounced? I'm implementing high resolution DOA algorithms including ESPRIT for radar systems. The word ESPRIT is obviously an acronym but, even so, how is it pronounced? Is the second i long like "ESPREEE" as the French would pronounce it or is it a short second i?

ESPRIT, like the French word.

For two reasons:

Out of vanity of myself, I learned it like that, I used it in presentations like that, and all I need for that to be correct is for a majority to recognize it to be correct. Hence, by stating this publicly, I make it correct! (just kidding, but when you look at language, there's no right or wrong – there's only conventions)
A. Paulraj himself pronounces it French in this presentation, and he should know :)

So, I disagree with @Envidia on how it's pronounced, but I agree with them on: It doesn't really matter.


How to partition collagen into separate fibers? I have a picture representing collagen fibers of a tissue:

Image is binary consisting of a background (pixels with value 0) and white collagen (pixels with value 1). A zoomed in section of an image looks like this:

I want to find a way to partition the white collagen into multiple fibers,
so that I could measure the length and width of each fiber separately.

One idea was to perform a
skeletonization and then split the skeleton on every branch point. The problem with this is that the irregular shape of the the original collagen induces unnecessary branches in the skeleton:

Here we see a skeleton composed of hundreds of segments, whereas the actual number of segments covering the "triangle" should be 3.

Any idea of an algorithm that could partition such image into multiple parts?

Your data seems to be quite noisy in a sense that there are a lot of holes in the fibers. When you know a typical minimal size of a fiber (probably its widths), you could use a binary closing filter:

Dilate the image with a structural element of the typical width of a fiber. Most likely a circle with a diameter equal to the minimal fiber widths could be appropriate.
Erode the image with the same structural element.
Use a skeletization algorithm on this data set.


xcorr in MATLAB for periodic function I have a periodic signal and I want to find it's autocorrelation function.
I can calculate it exactly:

$$R_{uu}(h) = \frac 1M \sum_{k=0}^{M-1} u(k)\cdot u(k-h)$$

But will xcorr() use this function for periodic signals?

BTW, I've created my own function for periodic signals:

function [R h] = intcor(u,y)
%Calculates the correlation between two vectors u, y
M=length(u);
h=0:1:M-1;
R=zeros(M,1);
for i=1:M
   for k=1:M
     if k-h(i) <= 0
     R(i)=R(i)+u(k)*y(k-h(i)+M); %since periodic, add M -> gets same result
   else
     R(i)=R(i)+u(k)*y(k-h(i));
   end
end
R(i)=R(i)/M;
end
end

The xcorr function assumes a linear cross-correlation, i.e. it assumes that the signals are zero outside the intervals. If you want to have an efficient implementation of the periodic cross-correlation, you can refer to the properties of the Fourier transform with

$$
\mathcal{F}(u \star y)=\mathcal{F}(u) \cdot \mathcal{F}(y)^*,
$$

while exploiting the fact, that the DFT considers signals to be periodic anyway. Here's the corresponding code for that:

N = 20;

u = randn(N,1);
y = randn(N,1);

R = intcor(u, y);

R2 = real(ifft(fft(u) .* conj(fft(y))))/(N);  % Calculate the periodic CC

hold off;
plot([u y R])
hold on;
plot(R2, 'ko');

legend({'u', 'y', 'intcor', 'fft-based'});

To answer your questions in the comments:

It is still not clear to me, if I can use xcorr to calculate the autocorrelation of a periodic signal correctly? -
No, you cannot use the xcorr function, as it assumes a linear cross-correlation. What you want is a calculation, where you assume the signal is periodic, and you just input one period of the signal.
What is the difference between linear cross-correlation and periodic cross-corellation? - Consider the general equation for cross-correlation:
$$
R_uy(h)=\sum_{k=-\infty}^{\infty}u(k)^*y(k-h)
$$

Normally, the summation is done over all times (i.e from $-\infty$ to $\infty$). However, the signals you pass to the xcorr-function are time-limited (i.e. they consist only of $N$ non-zero samples). Now, in the linear case, the cross-correlation assumes that the input-sequences vanish for $k<0$ or $k>N-1$. In the circular-correlation (i.e. assuming periodic signals), the signals are not assumed to be zero for the index above, but instead we assume $u(k+iN)=u(k), \forall i\in\mathbb{Z}$. I.e., the sequences are assumed to be periodic with period N.


What does MATLAB display when plotting a sound signal? When  you load the speech file 'mtlb' in MATLAB (which is by the way the voice of a woman saying the word 'MATLAB'), and plot the data, you get the following plot:

My question is this: what does the y-axis represent? What is MATLAB actually plotting here? Thanks!

What you are observing is the digital representation of the voltage, which in fact represents the acoustic pressure.
Workflow would be something like:

Vibrating larynx is producing Acoustic Pressure [Pa]
Variations of that pressure are converted by the microphone into Voltage [V]
Voltage is being sampled and quantized by the ADC (Analog to Digital Converter) and represented as series of numbers.

Now knowing your microphone's sensitivity [mV/Pa] you can relate the voltage and acoustic pressure. Additionally if ADC's parameters are known (or measured) then you can directly relate sample values to Pascals.


What is special about the theorem that over GF(2), $f^2(X) = f(X^2)$? What is special about the theorem that over GF(2), $f^2(X) = f(X^2)$?  There is a some 10-line proof for the same in "page no.42, Shu Lin, Daniel J. Costello-Error Control Coding (2nd Edition)-Prentice Hall (2004)"  I think the fact is more general and can be proved straight away.  I feel that  $f^m(X) = f(X^n)$ because in GF(2), $A^m=A=A^n$ where A is any polynomial, variable or constant from/over GF(2)and m and n are any two positive integers.  Am I missing some important point here?

You are confusing polynomials, elements of $\text{GF}(2)[X]$, with coefficients, elements of $\text{GF}(2)$.  Admittedly, it's certainly the case that as functions $\text{GF}(2)\to\text{GF}(2)$, any polynomial (or any function at all) will be equal to its (pointwise) power for any exponent greater than 0 which no doubt adds to the confusion.  Still, as polynomials (roughly, finite lists of coefficients), they need not be equal.  For example, $X^2 \neq X$ in $\text{GF}(2)[X]$ simply by definition.  With fields of characteristic $0$, evaluation is injective as a function $F[X]\to F^F$, so conflating polynomials with the functions they represent is not as much of a problem.

So for $f(X) = 1+X$, $$f(X)^3 = (1+X)^3 = 1 + 3X + 3X^2 + X^3 = 1 + X + X^2 + X^3 \neq 1+X^3 = f(X^3)$$


How to calculate mean and Variance of Gaussian-distributed random noise? I have been given a Gaussian-distributed random noise n(t) which has an average power of 5 mW.
How can I calculate the mean and variance of noise n(t)? Any suggestions? Thanks!

If the values for each sample of $n(t)$ are given, then the mean and variance can be estimated using the equations for sample mean and sample variance, which is trivial so I assume the OP is only given that it is Gaussian-distributed with average power of 5 mW. From that alone, there is no way to know what the mean of the signal is since it hasn't been specified that the process is "zero-mean" or not. If the signal does have a mean value, it will contribute to the total power so this must be specified.
The variance is the portion of the power that is not included in the mean. Since units of power were provided, then the resistance would need to also be given in order to provide units for any non-zero mean (voltage or current) if there is a mean value.


FIR Filter Deconvolution Suppose we have a system whose impulse response h has length K and fed with an input x that has length N. Then it is known that the output y has length M = K + N -1. This shows us the convolution matrix which relates x and y has size M x N. By existence theorem, the problems in the format $A_{m x n}x = b$ has a solution if and only if $m \leq n$. Since $m > n$ always holds, it seems to me that there will be always no solution for FIR filter deconvolution problem. Is this correct becase in a somewhere, I read about the possibility of existence.

A few comments:

An overdetermined system (with more equations than unknowns) can have exact solutions.
An overdetermined system can have approximate solutions; for example, in the least-squares sense, where $\mathbf{x}=(H^TH)^{-1}H^T\mathbf{y}$ minimizes $||H\mathbf{x}-\mathbf{y}||$.
Deconvolution is usually performed in the frequency domain, where $X(f) = Y(f)/H(f)$.


What Do We Expect Likelihood Function to Be Used For? In my understanding, I let make a example.

For example,  you can see in the following picture.

Consequently, We want to find  ^sigma. and We have already known the observation data(which is random data such as C) and bin(0.3).

To know "the estimated standard deviation ^sigma", First of all,

Find Likelihood for each bins with observed data( for example, data of C).
Find Maximum Likelihood among previous Likelihood values(from 1).
the Maximum Likelihood is "the estimated standard deviation ^sigma".  Am I wrong?
In other words, the Maximum Likelihood (^sigma) represent all bins's Likelihood values.

Normally you start with a probabilistic model of your data (observations). That model includes a noise term which has some variance or standard deviation which can vary over time or can be fixed. Then you write up the joint pdf of the data, normally written something like this $p(\boldsymbol{x}|\theta)$, where $\theta$ is a vector that includes all the parameters in the model including the noise variance. Then in likelihood estimation the pdf is viewed as a function of $\theta$ (data is fixed). The parameter vector that maximizes this function is the most likely parameter vector in the sense that it maximizes the probability that the model has generated the data. What your plots seem to indicate is that the underlying estimator is biased which means that the expected estimate is not the true value.

The ML estimator is a function. It takes as input your data $\boldsymbol{x}$ and provides a value (the estimate) $\hat{\theta}$. Your data is random, it is governed by a pdf. Because you plug random data into your estimator your estimate will also be a random variable governed by some pdf. I think this is what is illustrated in those images. Those plots illustrate the pdf of the estimate. It is then of course interesting to know if the expected value (the mean) of the estimate then is equal to the true value.


Filtering random Signal My question is easy one actually.

First, I generate a random signal using randn() function of MATLAB like this:

Then, I design a FIR filter of order 200 of pass-band characteristics with the pass-band $[0.2\pi, 0.4\pi]$ using the MATLAB function fir2():

My questions are:

What am I supposed to see when I filter a random signal using BPF?
What change in the characteristics will be occur considering the theory?

The output signal will still be normally distributed, but its power spectrum, i.e. its frequency content, will obviously be different from the input signal. If $S_X(\omega)$ is the power spectrum of the input signal, which is approximately flat, then the power spectrum of the output signal is

$$S_Y(\omega)=|H(\omega)|^2S_X(\omega)$$

where $H(\omega)$ is the frequency response of the filter.


Find which signal was transferred in 32QAM I've been assigned to solve a problem for 32QAM signals. The thing is, there are little to none sources in my native language to get any decent information on the topic.
The problem is (in its original form):

For 32QAM to determine which signal was transmitted using the minimum Euclidean distance between the symbols criteria, given that I and Q in the transmitted signal took values {+-1, +-3, +-5} and the received symbol sequence: [I, Q] = [1.5, -1], [-3.5, 4], [0.1, 0.1].

I assume, that [I, Q] can be replaced with [x, y] and the set of values is missing an A. If there are any misleading parts, I can specify them with my teacher, so you're welcome to ask.
It would be perfect if someone could explain the solution a bit for my further understanding or attach any sources to read.

Here's a sample 32QAM constellation.  With the scaling you were given, there's dots at every possible combination of $\lbrace -5 ... +5\rbrace$ except for the four at $(\pm 5, \pm 5)$.  "Minimum Euclidean distance" means that you find the dot that's closest to your actual received value and say that's the one -- i.e., if you receive $(-2.5, 1.2)$, then the nearest dot to that is $(-3, 1)$ and that's the one you decide on.
Just do that for the three I/Q pairs you've been given.  For extra credit, write a Python script...


How should I remove obvious anomalies from spectral data? I'm working with raman spectra and it seems that our spectrometer produces some anomalies at very specific wave numbers and I'm looking to remove them from the dataset. Here's a picture of 27 spectra of the same substance and below it a picture of the variance at the individual wave numbers.

From looking at the variance graph it seems fairly obvious where the anomalies are (ca. at wave numbers 480, 800, 1500, 1650-1700 and 2150).

I ultimately want to use the data for machine learning. So far I've implemented an algorithm that looks where the variance exceeds a certain threshold and then removes outlier spectra at that point (it removes the entire spectrum) until the variance doesn't exceed the threshold at any point. The problem with this algorithm is that it removes around 60% of the spectra which is acceptable but not ideal.

I'm thinking about a different method where I delete the information at all the wave numbers where the variance exceeds a certain threshold (like cutting out a vertical strip from the first graph where the variance is too high). Since most of the anomalies are relatively "narrow" I think I would't lose much information.

I want to know if one of these methods is clearly better than the other and why. Also I would like to know if there are any common methods used for dealing with this kind of data.

This seems to be a textbook case for a narrow median filter. Your $S$ is very broadband compared to your $N$, so a median filter will not affect your signal, as it can be designed quite narrow. And by using median instead of averaging (i.e. lowpass) you minimize the error introduced by the anomalies. I would give it a try.


DTFT of sawtooth wave through DTFT of rect signal In a course i'm currently taking, the lecturer computed DTFT for the following signal:
$$r[n] = \begin{cases}
1& 0 \le n \le N\\
0& \mbox{otherwise}
\end{cases}
$$
For $N = 32$ i pictured $\frac{1}{N}\frac{\sin{\omega*N/2}}{\sin{\omega/2}}$:

The lecturer told us that you can compute DTFT of any signal on finite support by using this formula:

\begin{align*}
\bar{X}(e^{j\omega}) = \sum_{k=0}^{N-1}X[k]\bar{R}(e^{\omega - \frac{2\pi}{N} k})\\
\end{align*}

So, I tried to render DTFT of sawtooth (red) as the sum:

That's clearly wrong. What did I do wrong?

The discrete Fourier transform (DFT) coefficients of a sequence $x[n]$ of length $N$ is defined by

$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N}\tag{1}$$

and the inverse transform is

$$x[n]=\frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{j2\pi nk/N}\tag{2}$$

The discrete-time Fourier transform (DTFT) of a sequence $x[n]$ of length $N$ is given by

$$X(e^{j\omega})=\sum_{n=0}^{N-1}x[n]e^{-jn\omega}\tag{3}$$

Plugging $(2)$ into $(3)$ gives

$$\begin{align}X(e^{j\omega})&=\sum_{n=0}^{N-1}\frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{-jn(\omega-2\pi k/N)}\\&=\sum_{k=0}^{N-1}X[k]\frac{1}{N}\sum_{n=0}^{N-1}e^{-jn(\omega-2\pi k/N)}\tag{4}\end{align}$$

With the DTFT of a rectangular sequence $r[n]=1/N$ for $n\in [0,N-1]$ and zero otherwise

$$R(e^{j\omega})=\frac{1}{N}\sum_{n=0}^{N-1}e^{-jn\omega}=e^{-j(N-1)\omega/2}\frac{\sin(N\omega/2)}{N\sin(\omega/2)}\tag{5}$$

Eq. $(4)$ can be written as

$$\sum_{k=0}^{N-1}X[k]R\left(e^{j(\omega-2\pi k/N)}\right)\tag{6}$$

So $(6)$ is the correct version of the formula, and the DFT coefficients $X[k]$ are defined by $(1)$. If you use those formulas correctly it should work.


Why Does Image Deconvolution Still Work with Image without Sharp Edges? Why does image deconvolution still work with image without sharp edges? Take this image for example:

Think of it in the frequency domain.
Some of the data is highly damped.
Namely energy in High Frequency is damped with low values.

Yet, unless it is zero identically "There is hope...".

Basically, this is the idea behind deconvolution algorithms.
Restore the energy in frequencies before it was damped.

Practically, we are limited in DR.
Namely if it is damped really hard, even if it not zero, in order to restore it one would amplify it significantly which will also lead to noise.

Regarding the edges, well, assuming we have the right way to recover the enrgies edges will be restored.


Linear time invariance RLC Why is a combination of R, L, C linear time-invariant system? But not diode, or transistor?

I looked over V-I equation for R, L, C, but still don't know how to prove it.

actually, the diode and transistor are also time-invariant components.  the transistor is even reasonably linear in a range of operation between cutoff and saturation.  even, for very small displacements, a diode can be sorta linear, but a diode in normal use is clearly not linear.


DoA estimation using microphone array I am using an array of 4 microphones in a square formation. My goal is to estimate the direction (azimuth only) of incoming speech sound.

The array looks like this:

I have recordings from the array to test my code with.

To estimate the delays between two microphone signals, I use cross-correlation method. This step is working for me and I've tested it somewhat thoroughly.

I am currently using delay values between $x1$ and $x2$, $x2$ and $x3$, and $x3$ and $x4$.

From the delays in samples, I calculate the azimuth angles of incoming sound between two microphones each using the arcus sine and the distance between the microphones. I get angles between $\pm\frac{\pi}{2}$. I am somewhat certain this step is working correctly too.

The next thing I do is mirroring the angles of either the first ($x1,x2$) or the third ($x3,x4$) combination based on whether the second angle (from $x2,x3$) is positive or negative.

From this, I should have the directions of the incoming sound, relative to the direction of $x2,x3$ being the front (right side of array in the picture).

To get the result of my DoA estimation, i add 90 degrees or $\pi/2$ to the second and 180 degrees or $\pi$ to the third angle. Now, all angles should point somewhat into the same direction, in my opinion, but it is not working.

Where is my mistake? I would really appreciate if someone could point out my error(s).

The mistake could be in assuming you know which angle $a$ to use when you don't know if the source is in front of or behind the array yet. In other words, you need to solve the same problem for $a$ that you're solving for $b$.

For your approach to work, you need to make sure side $a$ is the closest side to the source. You could find the closest side by finding the pair of mics where the maximum arrival time of the pair is less than the minimum arrival time of the opposite pair. This beam points outward, and the beam on the opposite side points inward. Only one pair of mics will meet those conditions at a time.

Also note that your approach for choosing which angle to use for $b$ only works if you consider the whole array as a single point (the distance between the mics relative to the distance to the source is vanishingly small). Otherwise, if the source is between the middle of side $a$ and one of $a$'s mics on the same axis as side $a$, angle $b$ will point the wrong way.

Instead of doing it like that, you might look at solutions using more than two mics at once.

An approach recently suggested here uses three equidistant mics and takes the cross product of the vectors pointing from one mic to another, with vectors having an additional component for the time difference.

Here's a really simple solution based on that, using axis-aligned diagonal mic pairs. Set your mics up like this:

      4
    1   2
      3

The time delta between mics 1 and 2 is the X component of a vector towards the source, and the time delta between mics 3 and 4 is the Y component. Check it out here.


Butterworth low pass filter zeros location after bilinear transformation explanation I am studying in a text book the transformation of a continuous time Butterworth low pass filter into a discrete time filter by means of bilinear transformation:
$$
    s = \frac{2}{T_d}*\frac{1-z^{-1}}{1+z^{-1}}
$$
which imposes the following relationship between discrete and continuous frequency variables:
$$
     w = 2*\arctan(\frac{\Omega*T_d}{2})
$$

As can be clearly seen from this formulas a value of $ \Omega $ = inf causes w = pi and a value of $ \Omega $  = -inf causes w = -pi. However, both values of w correspond to z = -1 in the Z plane

In the example a 6th order continuous time Butterworth low pass filter is designed. Its magnitude squared function in continuous time is defined by:
$$
  |H_c(j\Omega)|^2 = \frac{1}{1 + (\frac{j\Omega}{j\Omega_c})^{2N}   }
$$
From this formula it can be seen that as $\Omega$ goes to infinity the magnitude function goes to zero. So I can understand that in a 6th order filter there are 6 zeros in infinity (in S plane), which in discrete time means there are 6 zeros in w = pi and hence 6 zeros in z = -1. But:

What about the behavior of the system in s = -inf?
Isn't the magnitude squared function symmetric ?

Wouldn't it also cause the magnitude function to go to zero and should generate zeros on w = -pi and z = -1

An $N^{th}$-order analog prototype system results after bilinear transform in an $N^{th}$-order discrete-time system. So the number of zeros and poles remains the same. All zeros at $|s|\rightarrow\infty$ map to $z=-1$. In your case there will be $6$ zeros at $z=-1$ because there are $6$ zeros at $|s|\rightarrow\infty$.

Of course, the magnitude (squared) response is symmetric, but the zeros at $|s|\rightarrow\infty$ cause the zeros in the magnitude at $\Omega\rightarrow\infty$ as well as at $\Omega\rightarrow -\infty$ (you can imagine those points meet at infinity). This becomes more obvious in the $z$-domain where moving towards positive frequencies means moving counterclockwise from $z=1$ (i.e., $\omega=0$), whereas moving towards negative frequencies means moving clockwise from $z=1$. In both cases you'll end up at $z=-1$.


Differential Equation and Transfer Function of Multiple Summer Integrator Block Diagram I've been working on this more than 5 hour,
I know how to solve feedback, summer, but was not able to solve this one,
I've tried to simplify it by choose A(i) to be the first output of first summer and B(i) to be the second output of second summer

A(t) = u(t) + int(B(t))
B(t) = int(A(t) + y(t)
y(t) = int(int(B(t))

When I'm substituting A(t), its getting worse, forever recurs
y(t) = int(int(int(A(t))+y(t)))
y(t) = int(int(int(u(t)+int(B(t))

I've tried to use simulink, (I'm a newbie) and  I got 1 / (s^3 - s)
and I don't even know how to output the differential equation in matlab

How do I derive the solution?

First, transform the variables into Laplace domain for dealing with algebraic rather than differential equations, which greatly simplifies the labor. And then properly re-route those two feedback branches to simplify the block diagram yet still having the same overall transfer function.
All three sub-blocks of the system are integrators whose transfer function is $H(s) = 1/s$. And your feedback paths have unity gain.
Your initial system diagram looks like this:

The transfer function is defined like: $$ H(s) = \frac{Y(s)}{U(s)} $$
In the first step, lets move the upper feedback path, which is added to the output of the first integrator, to the left adder node.
Note that the feedback signal was $Y(s)$ . When you move this feedback signal to the front of the first integrator block, you should modify the feedback signal properly so that the two configurations produce exactly the same signals. Therefore, label the new feedback signal as $Y(s) s$, (more generally $Y(s)/H_1(s)$, where $H_1(s)=1/s$). This results in the following block diagram:

Observe that the middle adder is left free, and can be removed:

Then consider the lower feedback path which samples its signal from the input of the last integrator whose output is $Y(s)$. Therefore, the lower feedback path signal is $sY(s)$. Hence we can label the lower feedback path with the signal $sY(s)$.
This way, we can (apparently!) remove the feedback lines, as we have found their resulting (effective) input signals, and obtain the following final (simplified) block diagram in which you can merge the three cascaded blocks.

Finally, this is a simple input-output relationship (with feedback). : $$ Y(s) = \frac{1}{s^3} \{ 2sY(s) + U(s)\} $$
from which you find the transfer function as :
$$Y(s) \{ 1 - \frac{2}{s^2} \} = \frac{U(s)}{s^3} $$
$$H(s) = \frac{Y(s)}{U(s)} = \frac{1}{s^3 - 2s} $$


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

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

It appears that your understanding of harmonics is not entirely correct yet. Your example of two sinusoids with frequencies $f_1$ and $f_2$ can actually be explained in terms of harmonics IF
$$f_1=n_1f_0\quad\text{and}\quad f_2=n_2f_0,\qquad n_1,n_2\in\mathbb{Z^+}\tag{1}$$
is satisfied for some positive $f_0$. Note that this means that the ratio $f_1/f_2$ is rational. In that case the resulting sum is periodic. Choosing the smallest possible integer values of $n_1$ and $n_2$ such that $(1)$ is satisfied will allow you to compute the fundamental frequency $f_0$ of the resulting waveform.
In the case of your example we have $f_1=100$ and $f_2=201$ (units don't matter here). So we get $n_1=100$, $n_2=201$, and $f_0=1$ (same units as $f_1$ and $f_2$). So the sum is periodic with fundamental frequency $f_0=1$, even though the Fourier series coefficient of the fundamental is zero! Hence, we have a periodic signal with two harmonics and no fundamental. Nevertheless, the (fundamental) frequency of that signal equals $f_0=1$ (check that!).
With the period $T=1/f_0=1$ we have
$$x(t+T)=\sin\big[2\pi\cdot 100(t+T)\big]+\sin\big[2\pi\cdot 201(t+T)\big]=x(t)\tag{2}$$
where the last equality follows from the fact that we add an integer multiple of $2\pi$ to the argument of $\sin(x)$, which of course has period $2\pi$. Note that this doesn't show that $T$ is the smallest possible period, i.e., that $f_0$ is the largest possible fundamental frequency, but I trust that you can easily verify that.


Biquad filter from Matlab to C using "tf2sos" This is a biquad implementation I use:

float my_biquad(float in, float *coeffs, int scalefactor, float *delay)
{
    float out = coeffs[2] * in + delay[0];
    out = out * (float)pow( 2, scalefactor );
    delay[0] = delay[1] + coeffs[1]*in - coeffs[4]*out;
    delay[1] = coeffs[0]*in - coeffs[3]*out;

    // etc...

    return out;
}

It cascades several biquads:

for(int j=0;j<num_biquads;j++)
{
    tmp = biquad(...);
}

Each biquad has its coefficients and a scalefactor. Scalefactor is integer.

On the other hand, I have a Matlab filter:

[B_coeffs, A_coeffs ] = butter(4, 100/(16000/2), 'high');

I need to implement the Matlab filter using the above biquad.

I tried using tf2sos which gave me two biquad sections:

[ss, gn] = tf2sos(B_coeffs, A_coeffs);

but feeding ss coefficients to biquad in C gave different results than Matlab version. Gain gn is left unused, and I don't even have scalefactors.

How do I go from Matlab to this particular version of biquad in C? How to create proper coeffs and scalefactors?

What I can see from the code is the following mapping between the array coeffs[] and the filter coefficients $a_i$ and $b_i$ (note that $a_0=1$ and is not stored in coeffs[]):

$b_0=$coeffs[2], $b_1=$coeffs[1], $b_2=$coeffs[0]

$a_1=$coeffs[4], $a_2=$coeffs[3]

For the time being choose the scaling equal to $1$ (scalefactor=0) to check if everything else works. The final scaling per section will depend on the order of the sections, and on the goal of scaling (minimization of the possibility of overflow, or minimizing quantization noise).


What is the type number of a discrete time system given $H(z)$? Given a continuous time impulse response $h(t)$, if I take the Laplace transform and count the no. of poles at origin, that gives the type number of the system. For e.g., $$H(s) = \frac{2}{s(s+2)}$$ Here, there is one pole at origin and hence it is type 1 system. This also implies that the number of integrators in the system is 1. I'm having difficulty in determining the system type when it comes to discrete time systems. For e.g., given the $Z$ transform of impulse response, $$H(z) = \frac1z+\frac1{z^6}$$ how to deduce the system type? Is it 6 because there are six poles at origin? Also how does it relate to no. of integrators in the system given that it is discrete time. Does it correspond to no. of accumulators in any way? I came across this while searching but it talks about the order of the input polynomial as the type of the system which makes the type number independent of the impulse response. Is this correct?

I couldn't find a definition of system type for discrete-time control systems, but in analogy with continuous-time systems you should count the number of poles at $z=1$, not at $z=0$, because in discrete time, $z=1$ is the DC point, just like $s=0$ in continuous time.

A discrete-time integrator (an accumulator), has a pole at $z=1$. Its transfer function is given by

$$H(z)=\frac{z}{z-1}\tag{1}$$


What Are the Alternatives to FFT for Computing High Resolution Tone Power Levels? I have a system where there's a transceiver which transmits tones on specific frequencies (about 260kHz) and a receiver which is supposed to recognize those tones.

The transmitted tones are of low frequencies and has a very narrow band. At the receiver side I sample the band with ADC with sampling rate about 600kHz. The interesting band is about 1kHz wide, but whole analysed band is about 300kHz. To recognize if the signal exists I'm using 65k FFT? It takes about 100ms to collect samples. My goal is to find high-resolution tone power level.

I heard about zoom-fft and chirp-z, but is there another without expanding sampling time?

EDIT 1

Maybe i should ask a question diffrently:
What method (algorithm) gives the highest resolution (Hz) when I have a constant number of samples, and want to see only a small part of whole band (e.g. 1kHz in 300kHz)?

Actually all those "Angle of Arrival" algorithms can be used for that.

You can try:

MUSIC.
Pisarenko Harmonic Decomposition.
Sensor Array (MVDR).


Text recognition for Android app I´d like to use some OCR algorithm (Image to Text) into an Android app (Java). I have previously used several algorithms but for a PC (Matlab) app. Which library or algorithm should I use? Computational cost and response time are important issues.

As far as i know there are no native opensource Java OCR SDKs. There are Java APIs which wrap calls for native interfaces, for example, for one of the most popular opensource OCR engines - Tesseract (http://groups.google.com/group/tesseract-ocr/) - there are some Java wrappers like tesjeract (http://code.google.com/p/tesjeract/) or Tess4J (http://tess4j.sf.net/). That could work for you, but it's rather hard to set up and will require developing image-preprocessing and font training on your side.

One more solution could be a cloud service. It requires end-user application to have the internet connection, but it's independent from your programming language choice and resources limitations (which is importatnt on mobile devices, OCR proccess consumes rather big amount of recources). Have a look at www.ocrsdk.com for android ocr, it's a cloud-based OCR SDK that let you upload an image through web API and returns you the OCRed data.

This Web API based OCR SDK is not free, which may not be suitable for you, but i still recommend you try it out (it has a free 90 days trial without any upfront charges) as its pricing is really affordable in comparison with enterprise solutions while it provides enterprise-level OCR accuracy which is way better than open source.

You may also find useful this Android codesample at github (you need to get API key at http://ocrsdk.com to use it). Disclamer: i work @ ABBYY.


Pathloss effect MATLAB code I want to implement the pathloss effect using MATLAB simulation.
I use the D2D path loss model:
$$148+40\log_{10}(d)$$
I did not find the right way to add the pathloss effect in the received signal, which is a vector of QAM randomly chosen elements.

A propagation loss expressed in decibels is equivalent to an attenuation in linear scale. If you have a loss of l_dB (in dB), that's equivalent to an attenuation of l = 10^(l_dB/20). So, if your transmitted signal is S, the received signal is R = S ./ l.

Note that the factor of 1/20 in the conversion from dB to linear scale is needed because we're scaling the amplitude (voltage) of the signal.


A system having impulse response $ h(t)=u(t) $ stable or not? I know that for a system to be BIBO stable its impulse response must be absolutely integrable and the impulse response $ h(t)= u(t)$ integrates to approach infinity (i guess)
I proceeded as$$ \int_{-\infty}^{\infty}|u(t)|dt $$
$$ =\int_{0}^{\infty}dt $$
$$ =t\Big|_0^{\infty} $$
$$ =\infty$$
So the system must be unstable right? I m unsure  Can anyone explain in terms of poles of the $H(s)$?

Your derivation is correct. The impulse response is not absolutely integrable, hence the system is not BIBO stable. From the corresponding transfer function $H(s)=1/s$, you can see that there is a single pole at the origin. Systems with single poles on the imaginary axis, like the integrator in your example, are also called marginally stable. These systems are not BIBO stable, transients do not decay, and certain inputs can cause unbounded outputs.


How to find $h[n]$ system response of this equation? $2y[n-2]-2y[n-3]-4y[n-4]=x[n]-10x[n-1]-4x[n-2] + 4x[n-3]$  is the system that I'm looking for the response.

I transformed this system via using Z transform:

$$\frac{Y(z)}{X(z)}=H(z)=\frac{z^4 - 10z^3 - 4z^2 + 4z}{2( z^2 -z -2 )}$$

$$H(z) = 0.5z^2 -4z + \frac{-5}{z-2} + \frac{-1/3}{z-2} + \frac{1/3}{z+1}$$

How to find ROC to inverse the Z transfom to $h[n]$ ?

To take inverse Z-transform just the expression of H(z) is not enough. The inverse could be left sided or a right sided sequence depending on the causality of the system. Also stability of the system can give a clue. For example if the system is stable, the ROC will include the unit circle. Based on these additional information you should now find out right sided or left sided sequence from your partial function expansion. Your system seems to be anti-causal. So you can now work out the math and ask more questions if you get stuck.


Denoising / thresholding via wavelets When I apply differerent thresholding, wavelet denosing functions to non stationary time series which has been detrended via Loess regression and demean it. I expect that when this processed series are submited to denoising / thresholding will result in a clean series with smaller values than the submited signal and not to values at any point which are above those of the signal, Is my thinking correct.

On another hand could it be thought as per those values in a procesed signal which lie above the processed signal as if the procesed signal would ought  to be above those levels instead than below those levels. Of course the functions that described the signal is an aproximation so fiting erros should have to be expected.

I expect that when this processed series are submited to denoising / thresholding will result in a clean series with smaller values than the submited signal

No, you can get values that are greater. For example, consider the Fourier series of a signal. The Fourier series basis functions are wavelets. If we approximate the signal with only a few of the Fourier series coefficients you will get values that are above the original signal.

E.g. from the (image from wikipedia)


Evaluate frequency reponse of cascaded filter with bitwidth truncation in matlab How can we evaluate the frequency response of a series of fixed-point filters (H1,H2,...) cascading together in Matlab, such as:
12 bits input -> H1 -> bit truncation to 12 bits -> H2 -> bit truncation to 12 bits -> 12-bit output

As bit truncation is not a linear operation, there's strictly speaking no frequency response, as the overall system isn't linear.
But: assuming somewhat benign input, I'd simply model the quantization noise as uncorrelated (uff, I know that's a big jump!)
linear system H1 -> + quantization noise -> linear system H2 -> + quantization noise
or equivalently
linear system H1 -> linear system H2 ->  + colored quantization noise -> + quantization noise
allowing you to combine H1 and H2 and then just add non-white noise of the appropriate variance.
Note that the variances of your two noise sources depend on the expected quantization error, which in itself depend on the bit-widening properties of H1, H2. You could math out and find that variance, or you could run these filters with enough test signal and simply measure the squared error between 12-bit truncated and untruncated output.
I'm not a great fixed-point filter designer myself, but I'm pretty convinced there's plenty of material out there helping you figure out optimal truncation.


When to consider double (64 bit) floating point for Audio When synthesizing and processing audio on modern processors, when would one consider using something other than single precision (32 bit) floating point? Obviously, the audio coming in and going out to the real world is 16/24 bit, so I'm just talking about the precision of the signals (both the audio itself and things like filter coefficients) in the software.

Assume that:

the CPU/DSP has hardware floating point support for both single and double precision
The priority is quality audio, not high performance. For example, double precision would be considered if it offered better (perceptual) quality.

IEEE float singles only provide about 24 bits of mantissa.  But many DSP/filtering algorithms (IIR biquads with poles/zeros near the unit circle, etc.) require far more than 24 bits of mantissa for intermediate computational products (accumulators, etc.), just to get final results accurate to near 16 or 24 bits.  For these types of algorithms, 32, 40 and 48-bit scaled integer accumulators were often used with DSPs that had no FPU.

But on many current processor implementations (for PCs, smartphones, etc.), the double precision FPU is far faster than trying to use 32 or 64 bit scaled integer when your algorithm needs to have more than 24 bits of intermediate product.

To prevent trashing the data cache, the raw data can be in short integer or single precision float format, while only the more local computational kernel might use a higher resolution format.  But if you are sharing intermediate computation results between DSP modules, the interchange protocol between modules may also benefit from a higher resolution (more than 24-bit mantissa) bus or data format.


DFT N-point circular convolution of $\cos(2 \pi n/N)$ and $\sin(2 \pi n/N)$ I am attempting to solve the below problem:

I was able to calculate the N point DFT for a cosine which is similar to this (note, this is the DFT for $\cos(2\pi n k/N)$. So for this second image assume k in the image = 1 and m in the image =k).

I used a similar method to find the $N$ point DFT of $\sin(2\pi n/N)$ with the below results (same notes on the variables as above):

At this point I have $X_1(k)$ and $x_2(k)$, when I multiply these together I get a very ugly series that I can't solve. I think that I am missing something here, how should I solve this problem?

Thanks

HINT : $\cos(2\pi n/N) = 0.5e^{j2\pi n/N}+0.5e^{-j2\pi n/N}$. This itself is a IDFT equation, where $X_1[k] = 0.5\delta[k-1] + 0.5\delta[k+1]$.

Similarly for $\sin(2\pi n/N)$, $X_2[k] = (0.5/j)\delta[k-1]-(0.5/j)\delta[k+1]$.

So $Y[k] = (0.25/j)\delta[k-1] + (0.25/j)\delta[k+1]$, and $y[n] = (1/2) \times (0.5/j)e^{j2\pi n/N} -(1/2) \times (0.5/j)e^{-j2\pi n/N} = (1/2) \sin(2\pi n/N)$

To expand mathematically
$$
X_1[k] = \sum (0.5e^{j2\pi n/N}+0.5e^{-j2\pi n/N})e^{-j2\pi kn/N}$$ $$= \sum (0.5e^{j2\pi (1-k)n/N }) +(0.5e^{j2\pi (-1-k)n/N })$$ $$= 0.5\frac{(e^{j2\pi (1-k)N/N} - 1)}{(e^{j2\pi (1-k)/N} - 1)} + 0.5\frac{(e^{j2\pi (-1-k)N/N} - 1)}{(e^{j2\pi (-1-k)/N} - 1)}$$ $$= 0.5\frac{e^{j\pi(1-k)}\sin(\pi(1-k))}{(e^{j\pi/N(1-k)}\sin(\pi/N(1-k))} + 0.5\frac{e^{j\pi(-1-k)}\sin(\pi(-1-k))}{(e^{j\pi/N(-1-k)}\sin(\pi/N(-1-k))}
$$
You can use similar steps to compute $X_2[k]$. You need NOT proceed to multiply from here. Look at the numerator term of each fraction. $\sin(\pi(m-k))/\sin(\pi/N(m-k))$ is zero when $k \ne m$. For $k=m$, even before you reach the geometric series sum, you can see all terms in summation are 1. So $X[k] = N$ when $k=m$. So $X_1[k] = N\times 0.5\delta[k-1] + N \times 0.5\delta[k+1]$ . Though I am bit confused about the $N$ factor.


what's gaussian noise look like? I am using Matlab function imnoise to add gaussian noise to one image. However, it seems that adjacent samples of the noise is correlated. I am using matlab code below to add noise.

 I = imread('eight.tif');
 J = imnoise(I, 'gaussian', 0, 0.02);
 figure; imshow(I);
 figure; imshow(J);

From the image, it contains many short vertical gray lines which maybe be image structures for some cases. Is the noise really additive white gaussian noise(AWGN)?

histogram of noise:

Yes it is. However, unlike real AWGN which comes mostly from thermal agitation in electronics, matlab's noise is not purely random, this might be the cause of you seeing a deterministic pattern. But this pseudo-random noise is usually random enough to approximate real life noise

If you want to assess the noise "randomness", plot an intensity histogram of In-I, with In being the noised image and I the original image (you are subtracting the deterministic part and keeping the random part).


Why do Octave and MATLAB both automatically shift the gray values in this image? Edit: I found this other question that actually addressed a separate but possibly related problem I was experiencing. However, in the below case the image is in uint8 when read into both Octave and MATLAB so it doesn't actually address this specific problem.

I'm going through the Scientist and Engineer's Guide to Digital Signal Processing and experimenting with the image manipulation examples given later in the chapter. But I've run into an odd issue that is grinding my progress to a halt because it makes it impossible to predict what will happen.

I've used Octave from the start, but today downloaded a trial of MATLAB and it does the exact same thing.

Here's an example. I created a test image to experiment with grayscale transforms. Here is the original PNG:

In Octave I run the following code:

>> dots3 = imread('dots_3gray.png');
>> imshow(dots3);

And I get this as the result:

As you can see the values are all shifted -- the image as displayed by Octave (and also MATLAB) has downshifted all the values so the whites are now dark gray and the grays are now black. (ignore the size difference, the second is a screen capture) Both are 100x100 pixel grayscale images. The min value in the matrix as it is read in by Octave is 0 and the max value is 158.

Interestingly, using the eyedropper tool in Paint.NET I checked the range of values in the original PNG and the min value (in the center of a gray dot) is 96 and the max value (in the white field) is 255 -- a difference of 159 which is oddly enough exactly the number of steps in the range of numbers [0, 158] reported in dots3 as read by Octave. So some sort of shifting certainly appears to be going on here.

I can get the image to display with the correct colors by shifting the pixel values by 97, the difference between 255 and 158:

>> imshow(dots3 + 97)

Any idea why this behavior is happening? It makes it extremely difficult for me to experiment with the algorithms because I can't know if a change is caused by the code or by Octave / MATLAB doing this on its own.

Platform is Windows 10. Thanks.

Your image is an indexed image. Meaning it takes integer values that are supposed to be mapped to colours on screen via a particular colourmap.

You can see this if you read your image in as:

[I, M] = imread('dots_3gray.png');

where I becomes your 'indexed' image (which seems to contain indices from 0 to 158), and M is the colormap used to interpret it.

If you look at the documentation of imshow you will see that for indexed images specifically, you can specify which colourmap to use; in this case M to view it as expected:

imshow(I, M)

If you don't want this to be an indexed image, and you'd prefer to convert it to a standard grayscale image (i.e. in the range [0,1]) so that you can process it further, you can perform this conversion using the ind2gray command:

G = ind2gray (I,M); % convert to uint8 grayscale image in the range [0,255]
G = im2double (G);  % convert to double precision image in the range [0,1]

See my answer to a similar question in stackoverflow for more info.


Why FFT does not retrieve original amplitude when increasing signal length I don't understand why FFT return different maximum amplitude as the signal length increase. I would except that with a large signal length according to the frequency, detected amplitude will be very accurate.

From where this periodic signal is coming ?

Here is the related python code I used to generate the plot:

%matplotlib qt
%load_ext autoreload
%autoreload 2

import numpy as np
import matplotlib.pyplot as plt

def get_fft(x, dt):
    n = len(x)
    fft_output = np.fft.rfft(x)
    rfreqs = np.fft.rfftfreq(n, d=dt)
    fft_mag = [np.sqrt(i.real ** 2 + i.imag ** 2) / n for i in fft_output]
    return np.array(fft_mag), np.array(rfreqs)

def build_signal(amp, freq, signal_len):

    f = freq
    A = amp
    dt = 1
    t = signal_len

    x = np.arange(0, t, dt)
    y = A * np.cos(2*np.pi*f*x)

    fft_mag, rfreqs = get_fft(y, dt)
    return x, y, fft_mag, rfreqs, amp, freq, signal_len

# Build sin wave from 1 to 5000 signal length with freq=1e-3Hz and amplitude=0.5
all_t = []
all_amp = []
for t in np.arange(1, 5000, 100):
    x, y, fft_mag, rfreqs, amp, freq, signal_len = build_signal(amp=0.5, freq=0.001, signal_len=t)

    all_t.append(t)
    all_amp.append(fft_mag.max())

# Plot amplitude from FFT against signal length
plt.figure()
plt.plot(all_t, all_amp, 'o-')
plt.xlabel("Signal length")
plt.ylabel("Amplitude from FFT (0.25 is excepted)")

As MackTuesday ask in comment, I plotted the signal and related FFT for signal length = 1000 and 1400.

Your code is bit unclear, especially generation of your signal. Python allows for vectorized operations so it is good to use it. What's more, it is good to clearly specify the sampling frequency of your signal and use it then. Also please remember to normalize your FFT by length of your signal (in this particular case) and multiply by 2 (half of spectrum is removed so energy must be preserved).

Phenomena you are facing is connected with fact that your frequency is not matching exactly frequency bin for it. Because of that energy is leaking into other frequency bins. This is due to fact that you don't have integer amount of cycles of your sinusoid. Please search for spectral leakage. Easy way to check that is to use the equation for frequency spacing in frequency domain, that is:

$$\Delta f=\dfrac{f_s}{N} $$

Where $f_s $ is your sampling frequency. Then your frequency vector is:

$$f=[0, \; \Delta f, \; 2\Delta f, \; ...  \;,(N-2)\Delta f, \;(N-1)\Delta f] $$

If frequency you are analysing matches exactly one in your frequency vector, then amplitude will be correctly estimated, otherwise - leakage. In order to minimize this effect, you might want to apply window to your signal. It will reduce amount of energy leaking to other bins and in slight way will improve estimation of original amplitude, although you won't be really able to retrieve it in 100%.

Here is small modification of your code to help you visualize that. The red curve is ideal spectrum and blue one is for given number of frequency bins (strictly connected to your signal length).

import numpy as np
import matplotlib.pyplot as plt
import time

def get_fft(y, fs):
    """ Get the FFT of a given signal and corresponding frequency bins.

    Parameters:
        y  - signal
        fs - sampling frequency
    Returns:
        (mag, freq) - tuple of spectrum magitude and corresponding frequencies
    """
    n  = len(y)      # Get the signal length
    dt = 1/float(fs) # Get time resolution

    fft_output = np.fft.rfft(y)     # Perform real fft
    rfreqs = np.fft.rfftfreq(n, dt) # Calculatel frequency bins
    fft_mag = np.abs(fft_output)    # Take only magnitude of spectrum

    # Normalize the amplitude by number of bins and multiply by 2
    # because we removed second half of spectrum above the Nyqist frequency
    # and energy must be preserved
    fft_mag = fft_mag * 2 / n

    return np.array(fft_mag), np.array(rfreqs)

def build_signal(A, f, N, fs):
    """ Generate cosinusoidal signal of a given frequency and amplitude.

    Parameters:
        A   - amplitude of a signal
        f   - fundamental frequency
        N   - number of signal samples
        fs  - sampling frequency
    Returns:
        t   - time vector
        y   - generated signal
    """

    dt = 1/float(fs)            # Get the time domain spacing
    T  = N*dt               # Upper time limit
    t  = np.arange(0, T, dt)    # Construct time vector

    y = A * np.cos(2*np.pi*f*t) # Generate the signal

    return t, y

if __name__ == "__main__":

    # Define the parameters
    fs = 100   # Sampling frequency (must satisfy the Nyquist theorem)
    f  = 10    # Signal fundamental frequency
    A  = 0.5   # Signal amplitude

    # Prepare figure for plotting
    plt.figure()
    plt.ion()
    plt.show(block=False)

    f_off = 1                     # Frequency offset in Hz for zoom

    # Precomute the spectrum without any leakage
    (t, y) = build_signal(A, f, 1000, fs)
    (MAG, FREQ) = get_fft(y, fs)

    for N in np.arange(1000, 10010, 1):

        # Clear and scale
        plt.clf()

        # Generate the signal
        (t, y) = build_signal(A, f, N, fs)

        # Get the spectrum and frequencies
        (mag, freq) = get_fft(y, fs)

        # Plot the reference spectrum
        plt.plot(FREQ, MAG, '-or')
        plt.hold('on')
        # Plot the spectrum for given number of bins
        plt.plot(freq, mag, '-o')
        plt.title('Frequency spectrum for %d bins' % N)
        plt.xlabel('Frequency')
        plt.ylabel('Amplitude')
        plt.grid('on')
        plt.hold('off')
        plt.legend(['Ideal spectrum', 'Spectrum for %d bins' % N])
        plt.xlim((f-f_off, f+f_off))
        plt.ylim((0, A))

        plt.draw()
        plt.pause(0.5)

Two intermediate steps. Signal is 10 Hz sinusoid with an amplitude of 0.5, sampling frequency is 100 Hz.

Exact match:

Frequency component you are analysis falls between two frequency bins:


What factors affect the propagation delay at wireless communication? Generally, the signal propagation delay is d/c.
d is a distance between a transmitting antenna and a receiving antenna.
c is light speed.

However in NS-3, there is Random Propagation Delay Model

Random Propagation Delay Model: The propagation delay is completely random, and changes every time the model is called. All packages, even those that are sent between two fixed nodes, experience a random delay. As a consequence, the order of the transmitted packets is not maintained.

What kind of communication environment does this model take into account?

I think the obstacle, etc is considered. But because I can't find any official document or paper, I'm not sure about this.

Even if transmitter and receiver are fixed, you have unknown reflecting objects, resulting in unknown multi-path propagation. The resulting propagation delay can be modeled as a random variable. If there are moving objects in the environment (such as people moving in a room), then the propagation is not only unknown but also time-varying. The delay spread of a channel describes the difference between the maximum and the minimum propagation delay in the given environment.


Pywavelets CWT returning 0 after scale 64 I am running continuous wavelet transform in mexican hat wavelet basis on a data set (with millions of samples) using pywavelets. The command used is pywt.cwt(data, scales, wavelet). It works well for scales less than or equal to 64 but above 64 I am getting 0 as output for all coefficients. Any idea what is happening (I am assuming some truncation to 0 due to limited number of bits to store the data)? And how to fix it, if there is any solution.
Also above scale 64, the computation is getting really infeasible in terms of memory and runtimes, so even if I could fix it, I would actually prefer staying within scale 64. And of course in reality we would need infinite number of scaled daughter wavelets to complete the entire spectral range which is impossible. What are my other solutions, is it possible to have something like a low pass for the lower frequencies not covered using the mother and daughter wavelets, so I can visualize the continuous wavelet transform of high frequencies along with the short time Fourier transform of lower frequencies in one graph?

PyWavelets' cwt is flawed, and so is scipy's; use ssqueezepy.cwt.

a low pass for the lower frequencies not covered using the mother and daughter wavelets
... with the short time Fourier transform of lower frequencies in one graph

This is exactly what kymatio's Scattering1D does (pass average=False, oversampling=999), it's CWT for most frequencies, and STFT for lowest,
and one row dedicated to lowpass. Kymatio's CWT is identical to ssqueezepy's but with different filter design (more efficient for ML).


MACE correlation filter - any suggested values of of u? I am trying to understand the MACE correlation filter as per the paper

  Mahalanobis, A., Kumar, B.V.K.V., Casasent, D.: Minimum average correlation
  energy ﬁlters. Appl. Opt. 26, 3630–3633 (1987)

The equation to train the MACE filter is this

I however am abit confused about what values you use for u.  From my reading u, is an N x 1 vector containing the desired peak values for the training images.  Is this something I would just set to 1, could someone offer some suggested values

I have found a basic implementation in matlab from here (code posted below) but it does not give any definition to how I would go about deciding what u is?

D = diag(mean(abs(X),2));
% inv(A) * B = A \ B
XDX = ctranspose(X) * (D \ X);
h = (D \ X) * (XDX \ u);
H = reshape(h, size(I));

I will appreciate any help with this.

Thanks

The nth entry in the u vector corresponds to the value of the desired correlation peak of the nth training image. Typically you set u to be all ones if you have only images from the positive class. Entries corresponding to negative class example images are typically set to zero.


Understanding time delay estimation for GPS? I really need help to understand these questions which is highlight, for the GPS situation, I hope anyone who has been working in this situation before, please help me to have a clear understanding about "time delay estimation":

Here is the scenario:

Suppose I know precisely I have transmit the signal (data), I know the sampling, I know the time correspond to the very 1st point, I can go through the sign the time every single measurement, that I just transmitted. Now, I received the signal (data), I have section of data and in there somewhere is the signal that transmitted. Unlike, as the transmit, I know the time correspond to the very 1st point in my receive data:

what would I do to get the time delay?

what steps do we follow?

In there, somewhere, I don’t know where exactly, but in there somewhere is transmitted the signal and I want to estimate time delay between the transmit and receive. So, I want to look at:

what you need to do that?

What math. operation do I need to use to compute that time delay?

If I say:"I want to compute xcorr between 2 vectors". So, this xcorr is twice sum of the length – 1:

what do I do when I have xcorr output?

I take 2 vectors, and I compute even longer vector use xcorr, so now I compute the longer vector but I want a number, I don’t want a vector, I want a time delay estimation:

what do I do with the xcorr solved like a MATLAB?

If I say xcorr(x,y), x has certain length and y has certain length, the result is the length of x plus length of y and minus 1. So, now I have a long vector:

what do I do then to determine the time delay?

what the signal look like and what xcorr look like?

what do I do when I have xcorr to determine the time delay?

6.I know that have xcorr with bunch of complex number but I need only one real number. The real number that I want is the time delay. I take input data, I take received signal, I compute xcorr, then I do something to help myself to get the time delay. Result to that is single real number. But these are general ideas, I need to answer those question.

I have read this link: Time delay estimation of oscilloscope signals using cross correlation . It's really help my understanding, but they way he answered very high level, I think I need to understand the basic first, then continue with this link.

Please help me for these questions.

Thank you very much, everyone.

For GPS, (simplifying for now by omitting corrections for ionsphere and orbit position and relativistic clock offsets), we determine the "Pseudo-range" to each satellite (SV), which will be the relative delay between all the received satellites we have correlated to relative to our local clock- using correlation as you described (delay each locally generated SV until you get a maximum peak and that delay value is your pseudo-range. Thus, as long as we have a pseudorange for at least 4 satellites, we can solve a matrix involving 4 equations with 4 unknowns: x position, y position, z position and time (which is the time error between our local clock and the clock of the satellite). The satellite's time keeping is highly accurate (using Atomic clocks) and the 4 equation solution allows us to use very inexpensive clocks in our local receivers.

Below is shown a simplified 2D explanation showing how unique x,y position is determined on a plane with 3 satellites and perfect knowledge of time. Note how with 2 satellites there are two solutions. The same thing occurs with 4 satellites and unknown time, it is just that one of the solutions is on earth and the other can be neglected as it would be out in outer space.  In practice we benefit in a lower error solution by averaging over all satellites in view (forming a least mean square solution to an over determined set of equations).

Below is a sample output of what the output of a sliding correlator looks like for a particular SV in an actual captured GPS signal (the signal itself is indistinguishable from noise since the GPS signal prior to correlation is approximately 20 dB below the noise floor). Here we see the result of the spreading code that repeats every 1 ms, and the specific location of the peak in time relative to any local time reference we choose is the "pseudo-range".  These peaks are what you would get out of the Xcorr command (shown is the absolute value) and the horizontal axis where a peak occurs is the delay you are interested in finding. What is important is that all you need is the relative delay between all the satellites in view and your local clock reference (so you can pick any reference point you want as long as all codes are done with the same clock).  What I have seen in implementation, and perhaps this would be more complicated than you need, is to have a code NCO in a code tracking loop that is used to generate the codes for each satellite, and the code tracking is done using three correlators for each SV in an early-prompt-late arrangement each offset by 1/2 a chip. Early minus Late correlator output would then be proportional to code error (once acquired) and drives the Code NCO to track the delay. In such a fashion the code generator is always kept in alignment with the received code (in contrast to the sliding correlator approach I show below), and the relative delay is determined by the relative position of each code.

Below is the received signal prior to correlation (for the same signal used to generate the correlation plot above), with the blue representing the captured waveform after I had down-converted it to baseband and aligned it with the spreading code shown in red (as evidenced by the correlation peak). To note that any resemblance to the code that you think you might be able to see in the blue waveform to the red waveform is purely coincidental. The blue waveform in this view at this scale is visually indistinguishable from the AWGN process that is dominating it.

This is how the transmitted signal looks at that satellite (just two states +1 and -1):

And this is how that same signal looks as captured at our receiver, after down converting to baseband, removing all carrier and code offsets, and sampling at the mid point of each symbol (to match what we know was +1/-1 at the satellite). This matches very well to an AWGN noise process within the sampling space described, with no evidence of the signal without performing the correlation with the known sequence.

To add to the simplified descriptions above, here is a link to a good reference on the complete PVT solution (position, velocity and time) using GPS: http://www.nbmg.unr.edu/staff/pdfs/Blewitt%20Basics%20of%20gps.pdf


Total Harmonic Distortion THD for the rectifier I have a circuit with rectifier and this cause the non-linearity and consequently harmonic. I show it on the oscilloscope and then collect the signal from the oscilloscope to the computer by using the Labview. The procedure is to find the Total Harmonic Distortion. Could anyone please guide me what I have to do in order to find the THD?
I know that the f(input),f(sampling), N(number of samples)and M(number of cycles) are very important.
Thanks in advance!

this depends a lot how "good" this estimate needs to be. Here is a "rough & dirty" method.

Put your scope in spectrum analyzer mode. Make sure you see all harmonics.
Look where the highest harmonic is and that's still visibly sticking out over the noise floor.
Quadruple that frequency and use this as your sample rate.
Sample your signal for a sizable portion. You want at least a 1000 cycles
Filter the signal with a steep notch filter at the line frequency
Throw away the first 50 cycles or so.
Divide the power of the filtered signal by the power of the unfiltered signal. That's your THD estimate


What does natural scene image data-set mean? I came across a paper that is trying to classify and detect object(s) present in what they call a natural scene data-set. They have images containing objects like cars, bikes and people.

Can anyone tell me how this is different from other data-sets that contain the images of cars/bikes/people?

The site with the data-sets I am interested in is: Natural Scene Database.
For comparison, the Caltech Computational Vision data-set also contains cars and bikes from the rear and side view (with less background).

I want to know what is so special about the natural image data-set compared to the standard data-sets containing various viewpoints of cars/bikes from these set of images? Why would this natural image data-set be a bigger callenge in object detection?

From one perspective, not much. However, biological perceptual systems have evolved in response to physical properties of natural environments. And such environments possess certain characteristics, such as Laplacian distribution shaped histograms. Therefore natural images are grouped and analyzed together.


Are the FFT coefficients symmetric in image processing? On page 11 of Fundamentals of Image Processing by Ian T. Young, Jan J. Gerbrands, Lucas J. van Vliet (pdf) the results of the Fourier transform are shown (Figs 4a and 4b) and it appears to me (please correct me if I am wrong), that the transform may exhibit redundant data in four quadrants. If this is so, would it not be possible (especially for highly symmetrical images) to take just the top left quadrant, flip it left and then take the top half of this result and flip it downwards. Would this not preserve much of the information when a reverse FFT is performed? If so, it would allow as 4x reduction of transmitted information before a compression algorithm such as JPEG algorithm was applied.

For real images, there is indeed a formal redundancy, termed Hermitian or conjugate-symmetric as detailed by @Fat32. This symmetry however is "modulated" by the complex expression of the Fourier coefficients. So the FFT requires half the number of coefficients, but twice the amounts, due to the real/imaginary or modulus/phase couples. All in all, the redundancy vanishes.

More precisely, an FFT on a $2n$ length signal provides you with $n+1$ non redundant coefficients: the first (DC) and the last (Nyquist) are always real, the $n-1$ often complex. So you have $n-1$ complex coefficients, plus one pair of reals, altogether $n$ "pairs" of coefficients.

In general, with linear transformations, you cannot expect a bijective, one-to-one mapping with less transformed coefficients than samples in the signal/image, unless, as detailed by @Fat32, you reduce the data space, but then again, the transformed space would naturally have the same dimension as that of the data.


Power and Energy spectral densities of a linear system Here is a schematic of my desired circuit.

I want to calculate Power and Energy that the load dissipates during 0< t<8 You can see my attempts:

Therefore

How can find load's energy in Frequency domain for desired time
domain interval (0 < t < 8) ?
Did I obtain power correctly in
    frequency domain?

The relationships, with $x$ the input voltage to the network, $y$ the output voltage from the network, $R$ the load resistance, $p(t)$ the instantaneous power dissipated at the load, and $E$ the energy dissipated at the load and $t_0=8$

$$Y(f)=H(f)X(f)\\
y(t)=\mathscr{F}^{-1}\{Y(f)\}\\
p(t)=\frac{y^2(t)}{R}\\E\bigg|_{t=0}^{t_0}=\int_0^{t_0}p(t)dt=\int_0^{t_0}\frac{y^2(t)}{R} dt$$

IF the cable or network is a lossless cable (!), the example turns trivial:
$$y(t)=x(t)\\
p(t)=\frac{x^2(t)}{R}\\
E\bigg|_{t=0}^{t_0}=\int_0^{t_0}\frac{x^2(t)}{R} dt$$


How to measure the quality (or difference in quality) of different images? I would like to measure the quality of a color image, possibly as a difference to another image, but if possible, on its own too.

This question comes from the fact that we have a vision application, and improve lighting conditions, calbirations, and similar surrounding conditions. We want to measure the improvement of the quality of the image, such as:

Amount of detail
Contrast
Other things that describe the quality of an image or the amount of information in an image

Do you have ideas for this? Maybe histogram comparisons, or an advanced distance measure on the histogram? I'd love to hear some good ideas!

EDIT:

Of course the results of the vision application are a measure of the quality of the images. But:

it is not intended as a measurement for the quality, it is a measurement of something else
there are more parameters that influence the results of the application (for example the image itself), therefore it is only a rough measure of the quality of the image

I could use human evaluation, but I really want an objective measure, especially because the human eyes are not objective (our mind plays tricks, the monitor we watch it on influences the image, etc.).

FINAL EDIT:

After a discussion about what quality in an image is, I'll leave this question as is. My own ideas (and of my coworkers) are the following:

The range of the RGB channels (is it using all the colors possible)
The mean of the Sobel amplitude (the amount of detail or a rough measure of focus, also see this question about detail)
The amount of pixels that are clipped or 0 (no information in these pixels)
Also the mean of the RGB channels and the YCbCr channels (after conversion), to be able to see if the average color has changed (may lead to investigating light from other sources, or after changing our own lighting)

By combining these statistics we can get an impression of the quality from a more objective point of view, instead of just looking at the images. It also leads to pointers where the changes might come from.

It all depends on what you mean by "quality".  If it is visual quality as perceived by people, then the best thing to do is to ask multiple people to compare pairs of images to see which ones they like better.

Since you have mentioned a vision application, then it must have some performance metrics for whatever it does, which would be its measure of quality.

If by quality you mean how blurry or noisy the image is, then you might compute the Fourier transform to see whether the image contains high frequency components to see how blurry it is, or compute its entropy to see how noisy it is.


How do I implement the delay of the Q-Phase Input in GMSK? I am on a project to write a function in Java to perform GMSK modulation. I have a good understanding of the GMSK modulation block diagram to an extent, but I have a challenge on what it means to delay the Q-phase input by 90°.

Do I simply multiply my Q-phase bit stream by a negative sine carrier or what? Please, what is the correct thing to do?

A system that shifts the input signal by 90 degrees is a Hilbert transformer. For a single sine wave, a 90 degree phase shift is simple (e.g., sine becomes cosine), but for a general signal you need a Hilbert transformer. Check the basics here.
It can be implemented using an FIR (finite impulse response) digital filter. Such a filter can be designed in Matlab/Octave.

EDIT: The above is an answer to the question "how do I perform a $90^{\circ}$ phase shift (for non-sinusoidal signals)?". However, what I didn't notice - and what was luckily pointed out by Jim Clay in a comment - is the fact that in the given diagram there is no need for a phase shifter. The first $90^{\circ}$-block simply symbolizes the splitting of the signal in an I-component and a Q-component. It is not a phase splitter and the $90^{\circ}$-symbol is actually quite misleading.


Difference Between Two Forms of Equations of Auto Regressive (AR) Model I found equation 1 for Autoregressive model in various books and articles but I also found equation 2 for AR model, I understand the physical meaning of the equation but two different equations confused me a lot.

Can anyone tell me what is the difference between them and where one should use them?? I want to use AR equation for forecasting so which one should i use and why??. Both equations should yield different results because of the negative sign (which should just flip the actual signal i think so).

Auto Regressive Model means the current output is a linear combination of previous outputs and driving noise:

$$ y \left[ n \right] = \sum_{k = 1}^{p} {a}_{k} y \left[ n - k \right] + v \left[ k \right] $$

As you can see, the current value $ y \left[ n \right] $ depends on $ m $ values before it.
The parameter $ p $ is the order of the model while $ v \left[ k \right] $ is the driving IID White Noise.

Usually if we want to estimate the parameters of the model we have 2 main cases:

The Model Order $ p $ Is Known
In this case we need to estimate $ p $ parameters $ \left\{ {a}_{k} \right\}_{k = 1}^{p} $ from a set of $ m \geq p $ samples $ \left\{ y \left[ n \right] \right\}_{n = 1}^{m} $. the classic way to do so is using the Yule Walker Normal Equations.
The Model Order $ p $ Is Unknown
If $ p $ is unknown the intuitive a approach is estimate it and then go back to (1). One could estimate is by running on a grid of options (Try model order 3, 4, 5, 6, ... and chose the one which works best in MSE) yet this might do overfit.  A classic approach to evaluate the grid is using Akaike iInformation Criterion (AIC).

In your equations:

The connection is between the output to a different signal which means it is not an AR model.

If you tell us where did you get those form we might be able to clear things up.

Welcome to DSP!


Modeling multipath distortion and shifting by fractions of a sample I am working on modelling some FM transmitted signals, and I'm trying to figure out how to model the effects of multipath channels.  The most obvious method is to delay the signal by some integral number of samples and then add it to itself (with some attenuation factor).  But, this is obviously not very realistic, since 1 sample's worth of time is a long time for a radio wave (at least at a reasonable sample rate).  It would be nice to be able to delay the signal by a fraction of a sample.  However, I don't know how to do this...  Can anyone enlighten me?  I am fairly new to the DSP field, I'm sure this isn't a particularly challenging task.

Also, is this actually a good method for modelling multipath effects?

One approach is to upsample the signal, then delay it, then downsample (or, if feasible, you can do the rest of the processing at the higher sampling rate).

This general approach to modeling multipath is correct. However, you need to make sure that your model is realistic and/or useful for your purpose. A few things to consider:

If you add delayed/attenuated replicas of the signal to itself, then you're assuming you have line-of-sight between the transmit and receive antennas. The effect is different when you don't have line-of-sight.
The number of paths is important too. A channel with many paths is easier to understand theoretically, because you can make certain assumptions about the overall channel gain (e.g. that it's a Gaussian random variable).
The maximum delay is important too. A very large delay introduces very destructive effects on the received signal. Many, very small delays can be modeled as a single gain factor.
In a static environment, the channel is also static. If the antennas and/or obstacles in the environment are in movement, then the channel is dynamic, in the sense that the number of paths, and the gain of each, vary with time.


Transforming Non-Uniformly Sampled Image onto Rectangular Grid I have an image and associated latitude and longitude matrices. So I know the latitude and longitude of every pixel. However, due to the camera orientation, the image is skewed and thus the latitude and longitude sampling is non-linear. See the images below. I want to resample this image on rectangular, evenly sampled latitude-longitude grid. I have coded up a way to do this, but it's computationally expensive, and I feel certain that there must be a better way. I've done quite a bit of searching, but I haven't found good hits. Perhaps I just don't know the lingo necessary for the right search terms.

The image below is the original image.

I want to get a resampled image that looks something like this:

I produced this image, by basically moving a square boolean mask across the original image and finding the mean of the values as shown below:

Ideally, I think that pixels with mostly empty space should have an opacity value as well, but let's table that for now. Here's the MATLAB code I wrote to do this:

%% Create Skewed Image - All of These Are Knowns
rows = 1:size(frame, 1);
cols = 1:size(frame, 2);

% Create non-linearly sampled grid
lat_grid = zeros(size(frame));
lon_grid = zeros(size(frame));
for r = 1:numel(rows)
    for c = 1:numel(cols)
        lat_grid(r,c) = 0 + .001*rows(r) + 1e-3*cols(c) + 1e-6 * rows(r).^2;
        lon_grid(r,c) = 0 + .001*cols(c) + 1e-3*rows(r) + 1e-6 * cols(c).^2;
    end
end

h = figure; surf(lon_grid, lat_grid, frame, 'edgecolor','none'), view([0,90]),
colormap(gray(256))
set(gca, 'ydir', 'reverse')
ylabel('Latitude (deg)')
xlabel('Longitude (deg)')
hold on
filename = 'testAnimated.gif';

%% Create image in target grid
target_lat_array = 0:0.1:(max(lat_grid(:))+0.1);
target_lon_array = 0:0.1:(max(lon_grid(:))+0.1);
[target_lon_grid, target_lat_grid] = meshgrid(target_lon_array, target_lat_array);
resampled_frame = zeros(size(target_lon_grid));
n = 0;
for r = 2:size(target_lat_grid, 1)
    for c = 2:size(target_lon_grid, 2)
        n = n + 1;
        mask = ((lon_grid > target_lon_grid(r, c-1)) & (lon_grid < target_lon_grid(r, c))) & ((lat_grid > target_lat_grid(r-1, c)) & (lat_grid < target_lat_grid(r, c)));
        ph = patch(...
        [target_lon_grid(r, c-1), target_lon_grid(r, c-1), target_lon_grid(r, c), target_lon_grid(r, c)],...
        [target_lat_grid(r-1, c), target_lat_grid(r, c), target_lat_grid(r, c), target_lat_grid(r-1, c)],...
        [1e4, 1e4, 1e4, 1e4], [1,0,0,0.5]);
        if ~any(mask(:))
            delete(ph)
            continue
        end

        % Write to the GIF File
        im = frame2im(getframe(h));
        [imind,cm] = rgb2ind(im,256);
        if n == 1
            imwrite(imind,cm,filename,'gif', 'Loopcount',inf);
        else
            imwrite(imind,cm,filename,'gif','WriteMode','append');
        end

        delete(ph)
        resampled_frame(r-1, c-1) = mean(frame(mask));
    end
end

figure, imagesc(target_lon_array, target_lat_array, resampled_frame),
colormap(gray(256)), axis image
ylabel('Latitude (deg)')
xlabel('Longitude (deg)')

Anybody have any ideas? Like I said, I just feel certain that there must be industry-standard ways to do this.

You are trying to achieve something called interpolation and this can be done in roughly one line using MATLAB.

Interpolation is well explained over the web so I won't detail the functioning in this answer (Wikipedia page is very complete and redirect to many interpolation methods, here is the MATLAB documentation for 2D interpolation).

In a nutshell you parse as input both old locations (your longitude and latitude) and the data (the color corresponding to each location) as well as new locations (the rectangular and linear grid you want) and interpolation returns the new values corresponding to these locations.

Old locations could be either linear or non-linear as well as for the output locations. Moreover You can either ask for fewer or more points as the input number.

You should read carefully the documentation to address potential drawbacks of interpolation and to take account of the different methods available.

Note: if you have scattered data instead of gridded data you should use this function instead of interp2.

Here is an example with a non-linear sampled sine on 50 points.

The red plot shows a linear resampling that keeps the original number of points (same resolution) but some information is lost near $\pi/2$ as the new sampling rate is lower than the original one at this location.

The green curve illustrates a linear resampling with a greater number of points (150).This time less information is lost but more points are extrapolated between original ones (near $0.5$ and $2.5$).

Note that it is a trivial example to illustrate the possibilities given by interp1. it does not shows properly why information can be lost or extrapolated.


Direct downconversion (IQ downconverter) with ASK modulation-- what kind of Balun do you need? I have a simple but confused question.
I am trying demodulate a 915MHz ASK modulated signal with direct downconversion. So what I am building up is:
Antenna --> SAW bandpass filter @915MHz --> LNA --> Direct IQ down-converter (max2021 datasheet LO = 915MHZ) --> balun? --> decoder(sdr)
The max2021 downconverter is a 4-ports output: I+,I-,Q+,Q-. To put the differential signal back to single-ended, I was taught this is a balun's job.
But when I trying to build this system the balun concept is confusing. If successfully down-converted and DC blocked  , the IQ outputs are just unmodulated differential square pulses-- so what kind of balun do I use and at what frequency? Do I use a digital buffer in this case instead of an analog balun?
ps: I am using direct downconversion for carrier interference cancellation(move the carrier interference to DC then DC block) and design simplicity; and IQ downconverter to avoid phase ambiguity.

You could make your intentions clearer about "using direct downconversion for carrier interference cancellation" to "move the carrier interference to DC then DC block" in your scheme of "[d]irect downconversion ... with ASK modulation", had you attached the circuit drawing.
But, for a prescribed application of MAX2021 as GSM/CDMA signaling demodulator, you need a power combiner -- in fact, three power combiners: two for differential components in the $I$ and $Q$ paths and one to combine the single-ended $I,Q$ signals into one baseband signal. The "power" determinant should not confuse you: the name is inherited from a nomenclature of microwave power devices (e.g., microwave power amplifiers, attenuators, etc.), but, when used in the context of modulators/demodulators, designates exactly the component you need to combine $I_{+}/Q_{+}, I_{-}/Q_{-}$ outputs into single-ended signals. Besides, microwaves101.com directly indicates a correspondence between baluns and power combiners in the section 180 degree combiners.
Also, the baseband output of MAX2021 is not "square pulses" with zero rise/fall times; after output filtering according to a circuit shown in Figure 4. Baseband Port Typical Filtering and DC Return Network, baseband outputs look more like an "analog" signal with a frequency range defined by the symbol data rate.
In Figure 5. Demodulator Combining Diagram (page 15), the MAX2021 datasheet of your reference directs you to use Mini-Circuits' 180° power combiner ZFSCJ-2-1 and 0° power combiner ZFSC-2-1W-S+. ZFSCJ-2-1 features wideband operations, 1MHz to 3GHz. If your symbol data rate requires frequency below 1MHz, there are available on the market lumped-element/resistive power dividers/combiners of other suppliers, with a lower frequency value down to 4KHz and even DC.
Read also the MiniCircuits blog article and the application note AN10-006.


DFT Frequency domain analysis and interpolation I have a 2 part question, one may be related to why I'm not understanding the other.

A while back, I remember some professor saying that for the $N$ point DFT frequency domain, the values for $k$ around $N/2$ (halfway point for the spectral domain - those around $\pi$) are considered high frequencies, and the end points $(0, N-1)$ (those around 0, $2\pi$) are considered low frequencies.  I am having trouble seeing this mathematically given the definition for computing the coefficients of the DFT.

$$X[k]=\sum_{n=0}^{N-1}\left(x_ne^{\frac{-i2\pi{}kn}{N}}\right)$$

To me, it seems that the frequencies grow as $k$ increases right up until $k=N-1$, where the coefficients repeat due to periodicity.  Perhaps I am mis-remembering and this may be an incorrect assumption?

This leads on to my second question.  I have seen a few examples floating around online, where to interpolate in the time domain, one may pad in the frequency domain (running the IFFT after of course).  The examples I have seen seem to decide to pad the coefficients around the middle, i.e. for $N=4$ and a desired sampling of $2N$, [1 2 3 4] become [1 2 0 0 0 0 3 4].  Is there a reason why the middle is targeted and not the end like [1 2 3 4 0 0 0 0]?

If my first question is correct, then this simply pads the higher frequencies as 0.  I am aware that these zeroes add nothing to the frequency information as $X[k]=0$ cancels any effect during the inverse DFT.  But wouldn't moving any of the coefficients change some frequency information, i.e. for [1 2 3 4] becoming [1 2 0 0 0 0 3 4]

$N=4, k=2, X[k]=3$, becomes $2N=8, k=6, X[k]=3$

$$e^{\frac{+i2\pi{}2n}{N}}\neq{}e^{\frac{+i2\pi{}6n}{2N}}$$

A while back, I remember some professor saying that for the N point DFT Frequency Domain, the values for k around N/2 (halfway point for the spectral domain - those around π) are considered high frequencies, and the end points (0, N-1) (those around 0, 2π) are considered low frequencies. I am having trouble seeing this mathematically given the definition for computing the coefficients of the DFT.

You can define the DFT however you like, as long as in effect the same frequencies are covered - remember, $f(x) = e^{jx}$ is $2\pi$-periodic, so that, for example, the sum over $k=0,\ldots, N-1$ is the same as the sum over $k=-\left\lfloor\frac N2\right\rfloor,\ldots,0,\ldots,+\left\lfloor\frac {N-1}2\right\rfloor$.

  The examples I have seen seem to decide to pad the coefficients around the middle, i.e. for N=4 and a desired sampling of 2N, [1 2 3 4] become [1 2 0 0 0 0 3 4]. Is there a reason why the middle is targeted and not the end like [1 2 3 4 0 0 0 0]?

Again, different definitions of the DFT lead to different forms for zero-padding here.


LMS convergence and step size $\mu$ I am running the LMS algorithm based on Haykin's Adaptive filter theory.
I aim to plot the cost function $\mathbf{J}$ and calculate $\mathbf{J}_{\tt min}$ and the simulation
excess mean square error $\mathbf{J}_{\tt x}$.
I have two questions:

Based on theory, LMS converges when $0 \lt \mu \lt \frac{2}{\lambda_{\tt max}}$.
But this is not the
case here since while I have calculated theoretical $\mu=0.0019$, the
algorithm converges when $\mu$ is about $0.28$!!

Is this the proper way to calculate the simulation
excess mean square error $\mathbf{J}_{\tt x}$?

Could you help me find the mistake, either in my code or in theory?
Thanks for your time guys!!
clear all
clc
close all

sysorder = 3;
ap1 = 1.2;
ap2 = 0.53;
h=[1; ap1; ap2];

N = 1000;
N1 = 60;

lmw2 = zeros(1,N1);
lmw3 = zeros(1,N1);
lmw4 = zeros(1,N1);

x=randn(N,1);
corr_xx=xcorr(x);
for i=0:2
for j=0:2
R_xx(i+1,j+1)= corr_xx(N+i-j);
end
end

d=conv(x,h);

[R_dd,lags]=xcorr(d,1);
corr_xd = xcorr(d,x);
for i=0:2
R_dx(i+1) = corr_xd(N+i);
end

mu_max = 2/eigs(R_xx,1);%__________________________________________________mu step
%mu_max = 2/trace(R_xx);
mu = 0.1*mu_max;
%mu = 0.28;
%__________________________________________________________________________

[Jx,lmw,lme]=lm_s(x,d,N,N1,mu);

lm_mse = lme(1,sysorder:N).^2;%____________________________________________LMS MSE
Jmin = min(lme(1,sysorder:N).^2);%_________________________________________Jmin
Jx = mean(Jx);%____________________________________________________________Jexcess

figure(1)
plot(h, 'ko');
hold on
plot(lmw(:,N1), 'r*');

figure(2)
plot([sysorder-1:N1-1],lm_mse(sysorder:N1),'-','color','r');grid on;

%__________________________________________________________________________

function [Jx,wf,e]=lm_s(x,d,N,N1,mu)%______________________________________LMS algorithm
sysorder = 3;
w = zeros ( sysorder, 1 );
wf = zeros(length(w),N);

for n = sysorder : N1
   u = x(n:-1:n-sysorder+1);
   y(n)= w' * u;
   e(n) = d(n) - y(n);
   w = w + mu * u * e(n) ;
   wf(:,n) = w;
end

for n =  N1+1 : N
   u = x(n:-1:n-sysorder+1) ;
   y(n) = w' * u ;
   e(n) = d(n) - y(n) ;
   Jx(n) = w'*(u*u')*w;
end
end

LMS pseudocode:

Intitialization:
w [0] = 0

Computation:
for n = 0, 1, 2, 3, . . .
1. y[n] = wT[n]x[n].
2. e[n] = d[n] - y[n].
3. w[n + 1] = w[n] + µe[n]x[n].
end
```

In that range it is guaranteed to converge. It doesn't mean it will necesseraly won't converge for higher values. If you want deeper understanding you can read about the step size in Convex Optimization context where there the step size related to the Lipschitz Constant of the function (Which matches the eigen value for Quadratic functions).

If you share the problem itself from the book, we'll be able to solve the problem.


Detecting peak value on realtime IMU data stream I'm using Qt C++, QCustomPlot to read and display realtime value from an IMU. This is how I set the realtimeDataSlot:

void Settings::realtimeDataSlot(double x_acceleration_g, double y_acceleration_g, double z_acceleration_g, double z_acceleration_gnew)
{
    static QTime time(QTime::currentTime());
    // calculate two new data points:
    double key = time.elapsed()/1000.0; // time elapsed since start of demo, in seconds
    static double lastPointKey = 0;
    if (key-lastPointKey > 0.02) // at most add point every 20 ms
    {
      // add data to lines:
        ui->customPlot->graph(0)->addData(key, x_acceleration_g); // X axis
        ui->customPlot->graph(1)->addData(key, y_acceleration_g); // Y axis
        ui->customPlot->graph(2)->addData(key, z_acceleration_g); // Z axis
        ui->customPlot->graph(3)->addData(key, z_acceleration_gnew);

      lastPointKey = key;
    }
    // make key axis range scroll with the data (at a constant range size of 8):
    ui->customPlot->xAxis->setRange(key, 8, Qt::AlignRight);
    ui->customPlot->replot();

    // calculate frames per second:
    static double lastFpsKey;
    static int frameCount;
    ++frameCount;
    if (key-lastFpsKey >2) // average fps over 2 seconds
    {
      ui->statusbar->showMessage(
            QString("%1 FPS, Total Data points: %2")
            .arg(frameCount/(key-lastFpsKey), 0, 'f', 0)
            .arg(ui->customPlot->graph(0)->data()->size()+ui->customPlot->graph(1)->data()->size())
            , 0);
      lastFpsKey = key;
      frameCount = 0;
    }
}

which shows me as follows:

As a next step, I need to detect the peaks in any axis, say for example in the above figure in the Y axis there are peak values which I need to detect and count.

I marked in the peaks by hand in the figure. I define peak as the figure that value (positive values) more than 0.25 g at high rate.

Can somebody show me a way to do this realtime?

In order to detect peaks, the simplest method will be to take the difference of adjacent samples and check for zero crossing. You may also want to smooth out the data before taking difference to prevent spurious peaks due to noise. An outline of the procedure here would be like this

Smooth out the data using Moving average filtering. It is a simple filtering where you take average of, say, 5 neighboring samples to arrive at sample value.
For the above filtered signal take the adjacent sample difference ($y[n] = x[n] - x[n-1]) and save the sample indices where difference hits zero (or near zero).
To differentiate positive peaks from negative peaks, also check for sign of sample values around these indices. Eliminate negative peaks, and remaining ones indicate positions where positive peaks occur.
You may not get exact zero crossing, so you may want to keep a threshold. Also, the position may have offset due to filter delay. A 5-sample MA filter will delay your sample train by ceil(5)=3 samples.
There are ways to combine MA filtering and taking difference. Look for filters which do differentiation (https://www.dsprelated.com/showarticle/35.php)


Measuring distance of periodic structure in biological sample with FFT in MATLAB I have pictures, taken with a microscope, of biological samples that exhibit a periodic pattern. Previously, I have used ImageJ to manually rotate the image perpendicular to the pattern and plot its grey value profile. Doing FFT in MATLAB I can determine the distance between two peaks of that pattern via its spatial frequency. But this all depends on how (well) I rotate the sample. So my goal is to actually use fft2 in MATLAB to directly transform the image regardless of sample angle. Somewhere I have read about calculating the rotational average (after fftshift) of the resulting spectrum, but I don't know how.

Is this the right approach anyway?
Or is there an easier way to calculate the (average) distance between two units of a periodic pattern in an image without having to adjust the angle beforehand?

i am not an image processing guy.  but i know what the 2-dimensional Discrete Fourier Transform is.

normally, in two dimensions, the way we measure period is with an operation related to the autocorrelation and you can get to the autocorrelation $R_x[k,\ell]$ with the DFT and magnitude-squaring the result $|X[m,n]|^2$ and inverse DFT to get the autocorrelation.  mathematically, you should have half of the 2-dim buffer (like quadrants II and IV) set to zero, to prevent overlap in the correlation, then you have a triangular envelope on the peak values, even if it were a perfectly periodic function.

so, in the $k$ dimension of the result there will be an apparent period, call it $p_k$ and similar for the $\ell$ dimension.  the period will be $p = \sqrt{p_k^2 + p_\ell^2}$ .  actually, when viewing $R_x[k,\ell]$ as a surface, there should be a circular ridge of radius $p$.


FFT and number of samples relations I am new in signal processing. I generated a signal with $f_{in}=10 \mathrm{kHz}$ and also take 64 sample from this signal after doing some process in an ADC block. I want to convert the result to the main signal and compare them whit each other. So if what is the difference between taking an FFT of the 64 points or other numbers like 1024 in the same measured time and what should be changed?

Thank you

It is not clear from your question why you are even taking an FFT; knowing that may also lead to additional insights- but the main point to understand is frequency resolution in the FFT; assuming you are not doing any windowing (shaping the envelope of the time domain signal before taking the FFT), the frequency resolution is 1/T where T is the length of your FFT block in seconds. From this you will see that the frequency resolution is one FFT bin, and is equal to $f_s/N$ where $f_s$ is the sampling rate used and N is the number of samples in the FFT.

For additional info on this see What happens when N increases in N-point DFT and the additional referenced posts there.