what is the relationship between the reliability of a monte carlo test and the length of the samples tested and the number of errors found?

For each sample, we note $1$ if error and $0$ otherwise. Then, for bit error rate $q$, the tests are independent Bernoulli random variables $\{X_i\}$ with probability $q$.
The estimate is $\hat{p} = \frac{S}{n}$  where
$$S = \sum^n_i X_i \tag{1}$$
and $\mu = \mathbb{E}[X] = np$.
Chernoff bound (see Corollary 5) tells us that for every $0 < t < \mu$,
$$\mathbb{P}\left( \mid S-np\mid \leq \frac{t}{\mu}\mu\right) \geq 1 - 2 e^{\frac{-t^2}{3\mu}} \tag{2}$$
Choose $t=\epsilon np$ with $0 < \epsilon < 1$,
$$\mathbb{P}\left( \mid S-np\mid \leq \epsilon np\right) \geq 1 - 2 e^{\frac{-\epsilon^2 np}{3}} \tag{3}$$
$$\mathbb{P}\left(p(1-\epsilon) \leq  S/n \leq p(1+\epsilon)\right) \geq 1 - 2 e^{\frac{-\epsilon^2 np}{3}} \tag{4}$$
Therefore, $\epsilon$ is relative estimation error and the right hand side of $(4)$ is confidence level which depends on sample size $n$.
For example, if you want estimate a BER within $10^{-3}\times(1 \pm 10\%)$, using sample size $n=8.99 \times 10^5$ will guarantee a confidence level greater than 90%.
Note that choosing another $t$ and/or using another bound may result in a tighter required sample size.


How to scale a signal to get desired variance

standardizing a Normal random variable $x$
$$
\frac{x -m}{\sigma} \sim \mathcal{N}(0,1)
$$
so if $y$ is a random variable that is Normal with mean zero, with standard deviation $1$,
$$
\sqrt{\sigma^2}y \sim \mathcal{N}(0, \sigma^2)
$$
This standardization actually works for more than just Normal Distributions.


How to determine range of HSV values of the image?

You can simply specify a ROI for that region and convert it into HSV.

For eg (below is pseudo-code in Python-OpenCV)

# define ROI of RGB image 'img'
roi = img[r1:r2, c1:c2]

# convert it into HSV
hsv = cv2.cvtColor(roi,cv2.COLOR_BGR2HSV)

Now it gives you the hsv values of the region.

But one or two difficulties there:

Your object may comprise some part of the ROI. There may be some background region also inside your ROI. ie if your object is red and background is blue, your ROI will contain two different values, one for RED and other for BLUE. Which one will your program select?
But colors won't be constant like RED or BLUE. You get wide variety of colours. Below is a plot of different colours for a ROI with green object.

But it has values from yellow to slight blue with green as major.

So you will have to select color with maximum probability, may be from checking its histogram.

Check this SOF for more details : https://stackoverflow.com/questions/13111504/selecting-best-range-of-values-from-histogram-curve


Why entropy is undefined in low contrast image

You are getting a nan returned because you are performing an operation that yields undefined result.

Based on your code, the only option I can see is the log operater that theoretically yields $-\infty$ when you attempt to take the logarithm of a 0-valued pixel. C++ will thus consider it not a number. Other functions I have seen for calculating entropy such as this example in octave, just ignore the zero pixels to solve these cases.


How do can we cancel out noises while listening voice on headphones?

here is a great website on building noise cancelling headphones http://headwize.com/?page_id=684

The concept behind noise cancelling headphones is simple. You listen to the environment and then play it back 180 degrees out of phase. When the original wave and the shifted wave meet, they cancel eachother out

The circuit is pretty simple. You have microphone go into an inverting amplifier, then you play it through the speaker (simplified version) This achieves noise cancellation.

edit

one major benefit of this solution is it is purely implemented in (analog) circuits. It doesn't require adaptive filtering, a decision engine, or any computing power. This type of circuit is what most commercial products use: Active filtering

Sometimes the solution to the problem isn't what we expected. We can get into signal estimation, adaptive filtering, DSP boards or MCUs, or any number of complex mathematically rigorous constructs but really, all that's needed is an inverting amplifier. If you want a software solution please state so (and in that even we can go back to signal estimation and adaptive filters) but since you asked for practical solutions to actually convert regular headphones I have presented a hardware implementation.

edit2

unfortunately the original website seems to have vanished, but here is the same exact content in a PDF I include the full article because it has parts lists, some helpful troubleshooting tips, and far more detail than I've included. And in case the new link breaks here is the main circuit with a brief explanation from pdf. The circuit itself looks complex but its really a 3 stage amplifier but one for each audio channel (left and right) so 6 amps total. As you can see the second stage is the inverting amp (makes signal 180 degrees out of phase) achieve our noise cancellation.

The following material is an excerpt from the noise cancellation headphone guide written by Jules Ryckebusch  and all credit goes to Jules for design and explanation.


Inverse Discrete Time Fourier Transform of $1$

The IDTFT of $X(e^{j\omega})=1$ is indeed

$$x[n]=\frac{\sin(n\pi)}{n\pi}\tag{1}$$

Now, what happens for indices $n\neq 0$? As it turns out, you can safely rewrite $(1)$ as

$$x[n]=\delta[n]\tag{2}$$

where $\delta[n]$ is the discrete-time unit impulse. (HINT: think about where the zeros of $\sin(x)$ are).


Counting shape changing pattern in signal

You can trying using a musical pitch detection/estimation algorithm, which may be more robust in some cases than just a simple autocorrelation (especially in one of your examples that includes a slightly inharmonic overtone).  From the frequency or period of the estimated "pitch" you can compute the number of pattern repetitions per second.

A bunch of pitch estimation methods are mentioned in my answer to this question: https://stackoverflow.com/questions/7181630/fft-on-iphone-to-ignore-background-noise-and-find-lower-pitches


How to do continuous signal processing (i.e without windowing)?

Yes, it's possible to analyse sound the way ears do.

For example, you could compute the DFT of a signal continuously using several Goertzel filters.

$$
y_k[n] = e^{j2\pi k/N} y_k[n-1] + x[n]
$$

where $k= 0,1,\ldots, N-1$, so that $y_0$ is the DC or zero frequency term.

Of course, this is an unstable filter, so some resetting or forgetting factor is needed to keep it stable. There are ways to make it better behaved numerically (see pp84-93 of Stoer & Bulirsch).

The reason why we don't generally do things this way is that we are usually trying to do them in the most computationally efficient way possible.  And nature doesn't necessarily do this, so we find ways to do similar things but that work better on a computer.


Gaussian noise generation for a given SNR ?

Since nobody answered this question yet, i'm going to do my best to do it. Be aware i'm no expert on DSP.

I believe your first function is calculating the signal and noise average power based on each one's Fourier Transforms. So you're returning the SNR in decibels at the end.

\begin{equation}
SNR_{db} = 10 * \log_{10}{\frac{P_{signal}}{P_{noise}}}
\end{equation}

Your gaussian noise function generates the noise based on a scaling factor k of the signal max amplitude. Since you want to scale the amplitude of the noise based on your signal, i believe you want a relationship of:

\begin{equation}
k=\frac{A_{noise}}{A_{signal}}
\end{equation}

With each A meaning RMS amplitude. That can be generated used the following equation:

\begin{equation}
k = \sqrt{\frac{1}{SNR_{linear}}}
\end{equation}

or, to be more clear:

\begin{equation}
k = \sqrt{\frac{P_{noise}}{P_{signal}}} =\sqrt{\frac{A_{noise}^2}{A_{signal}^2}}
\end{equation}

If you were given the SNR in decibels and was asked to generate a noise based on it, you can use the following equation:

\begin{equation}
k = \frac{1}{10^{\frac{SNR_{db}}{10}}}
\end{equation}

In the third function you're generating the output signal by adding the frequency components of each signal, but if it's just an additive gaussian noise, you could just add the noise to the signal. (I'm not exactly sure on this).

You could also generate the linear SNR from your SNR in decibels, I've used this function in one of my projects once:

def linear_snr_from_db(snrdb):
    return 10.0 ** (snrdb/10.0)

I got all of the equations based on the wikipedia page for SNR and my little bit of experience.


Fourier Transform within a certain Limit

So, your integral with bounds

$$X(f) = \int_{a}^b x(t) e^{-i2\pi ft}\,dt$$

could really be written as

$$X(f) = \int_{-\infty}^\infty (r_{a,b}(t)\cdot x(t)) e^{-i2\pi ft}\,dt$$

with $r_{a,b}(t) = \begin{cases}1 & a<t\le b\\0&\text{else,}\end{cases}$  i.e. as a "windowed" view at $x(t)$.

Luckily, the convolution theorem of the Fourier transform allows us to say that

$$\mathcal F\{r_{a,b}\cdot x\} = \mathcal F\{r_{a,b}\}*F\{x\}$$
with $*$ being the convolution operator.

Thus, you'd just need to numerically determine an estimate for the Fourier transform of $x$, $F\{x\}$, and convolve that with the Fourier transform of your window $r$. Luckily, the Fourier transform of a unity-wide rectangle is known (it's a sinc function), and you just need to apply the scale and shift theorems of the Fourier transform to convert that standard sinc to the Fourier transform $F\{r_{a,b}\}$.

Then you'd need to estimate the convolution of the two.

Assuming your signal is well-behaved (especially: band-limited), a discretized signal $x[n]$ represents your $x(t)$, and hence, you could estimate the continuous time Fourier Transform within the signal bandwidth by doing a DFT (matlab: fft()) of sufficient samples of a period of $x(t)$.

As said, finding $F\{r_{a,b}\}$ is easy to do analytically. Simply generate the same amount of samples of $F\{r_{a,b}\}$ as you did for $x(t)$, and numerically convolve the result (matlab: convolve, filter).

This might be much more effort or much less elegant or much less accurate than simply doing numeric integration, which matlab can also do out of the box. But I'm assuming you're not interested in that, because a) it's trivial to find "matlab numerical integration" if you just look for it, so I presume you already did that and b) I'm assuming you're asking how to solve this with signal processing methods, considering the place you ask.


For complex values, why use complex conjugate in convolution?

Turns out that convolution and correlation are closely related.  For real signals (and finite energy signals):

Convolution: $\qquad y[n] \triangleq h[n]*x[n] = \sum\limits_{m=-\infty}^{\infty} h[n-m] \, x[m]$

Correlation: $\qquad R_{yx}[n] \triangleq \sum\limits_{m=-\infty}^{\infty} y[n+m] \, x[m] = y[-n]*x[m]$

Now, in metric spaces, we like to use this notation:

$$ R_{xy}[n] \triangleq \Big\langle x[m], y[n+m] \Big\rangle = \sum\limits_{m=-\infty}^{\infty}x[m] y[n+m] $$

The $\langle \mathbf{x},\mathbf{y} \rangle$ is the Inner product of the vectors $\mathbf{x}$ and $\mathbf{y}$ where $\mathbf{x} =\{x[n]\}$ and $\mathbf{y} =\{y[n]\}$.  Then we also like to define the norm of a vector as

$$\begin{align}
 \| \mathbf{x} \| &\triangleq \sqrt{\big\langle \mathbf{x},\mathbf{x} \big\rangle} \\
 &= \sqrt{\sum\limits_{m=-\infty}^{\infty}x[m] x[m]} \\
 &= \sqrt{\sum\limits_{m=-\infty}^{\infty}x^2[m]} \\
\end{align}$$

and that looks a lot like the Euclidian length of a vector with an infinite number of dimensions.  All this works very well for the case where the elements $x[n]$ of the vector $\mathbf{x}$ are all real.  The norm $\| \mathbf{x} \|$ is always real and non-negative.

So, if we generalize and allow the elements of $\mathbf{x}$ to be complex-valued, then if the same definition of norm is to be used,

$$\| \mathbf{x} \| \triangleq \sqrt{\big\langle \mathbf{x},\mathbf{x} \big\rangle}$$

then the definition of the inner product needs to be modified a little:

$$\Big\langle \mathbf{x},\mathbf{y} \Big\rangle = \sum\limits_{m=-\infty}^{\infty}x[m] y^*[m]$$

Then if $\mathbf{x}$ has complex-valued elements, the norm comes out as:

$$\begin{align}
 \| \mathbf{x} \| &\triangleq \sqrt{\big\langle \mathbf{x},\mathbf{x} \big\rangle} \\
 &= \sqrt{\sum\limits_{m=-\infty}^{\infty}x[m] x^*[m]} \\
 &= \sqrt{\sum\limits_{m=-\infty}^{\infty}\Big|x[m]\Big|^2} \\
\end{align}$$

So, evidently, Haykin is just walking that definition of inner product back to the definition of convolution.


What is difference between repeated poles and distinct poles?

Distinct poles do not need to share an x or y coordinate. They are classified as distinct as long as they do not share the same x AND y space. So any two poles are either distinct or repeated.

The effect repeated poles have on the impulse response of a filter is a little complicated, but the short answer is that it does change. This brief article describes the effects of repeated poles as well as providing proofs:

https://www.dsprelated.com/freebooks/filters/Repeated_Poles.html

Hope that helped.


FFT equivalent for generalized unitary transforms

It's all about structure. One early paper on this is  A Unified Treatment of Discrete Fast Unitary Transforms, 1977:

  A set of recursive rules which generate unitary transforms with a
  fast algorithm (FUT) are presented. For each rule, simple relations
  give the number of elementary operations required by the fast
  algorithm. The common Fourier, Walsh-Hadamard (W-H), Haar, and Slant
  transforms are expressed with these rules. The framework developed
  allows the introduction of generalized transforms which include all
  common.transforms in a large class of "identical computation
  transforms". A systematic and unified view is provided for unitary
  transforms which have appeared in the literature. This approach leads
  to a number of new transforms of potential interest. Generalization to
  complex and multidimensional unitary transforms is considered and some
  structural relations between transforms are established.

Among the most common, the discrete sine (DST), cosine (DCT), Hartley, wavelet transforms and many others (Walsh, Hadamard, Paley or Waleymard, triangle, jacket, slant, Hermite) have fast counterparts.

One global initiative towards a systematic construction was termed Algebraic Signal Processing Theory. Let me mention two papers laying some foundation using basic algebraic structures:

Algebraic Signal Processing Theory: Foundation and 1-D Time, 2008

  This paper introduces a general and axiomatic approach to linear
  signal processing (SP) that we refer to as the algebraic signal
  processing theory (ASP). Basic to ASP is the linear signal model
  defined as a triple ($\mathcal{A}$, $\mathcal{M}$, $\Phi$) where
  familiar concepts like the filter space and the signal space are cast
  as an algebra $\mathcal{A}$ and a module $\mathcal{M}$, respectively.
  The mapping $\Phi$ generalizes the concept of $z$-transform to
  bijective linear mappings from a vector space of signal samples into
  the module $\mathcal{M}$. Common concepts like filtering, spectrum, or
  Fourier transform have their equivalent counterparts in ASP. Once
  these concepts and their properties are defined and understood in the
  context of ASP, they remain true and apply to specific instantiations
  of the ASP signal model. For example, to develop signal processing
  theories for infinite and finite discrete time signals, for infinite
  or finite discrete space signals, or for multidimensional signals, we
  need only to instantiate the ASP signal model to a signal model that
  makes sense for that specific class of signals. Filtering, spectrum,
  Fourier transform, and other notions follow then from the
  corresponding ASP concepts. Similarly, common assumptions in SP
  translate into requirements on the ASP signal model. For example,
  shift-invariance is equivalent to $\mathcal{A}$ being commutative. For
  finite (duration) signals shift invariance then restricts
  $\mathcal{A}$ to polynomial algebras. We explain how to design signal
  models from the specification of a special filter, the shift. The
  paper illustrates the general ASP theory with the standard time shift,
  presenting a unique signal model for infinite time and several signal
  models for finite time. The latter models illustrate the role played
  by boundary conditions and recover the discrete Fourier transform
  (DFT) and its variants as associated Fourier transforms. Finally, ASP
  provides a systematic methodology to derive fast algorithms for linear
  transforms. This topic and the application of ASP to space dependent
  signals and to multidimensional signals are pursued in companion
  papers.

Algebraic Signal Processing Theory: Cooley-Tukey Type Algorithms for DCTs and DSTs, 2008

  This paper presents a systematic methodology to derive and classify
  fast algorithms for linear transforms. The approach is based on the
  algebraic signal processing theory. This means that the algorithms are
  not derived by manipulating the entries of transform matrices, but by
  a stepwise decomposition of the associated signal models, or
  polynomial algebras. This decomposition is based on two generic
  methods or algebraic principles that generalize the well-known
  Cooley-Tukey FFT and make the algorithms' derivations concise and
  transparent. Application to the 16 discrete cosine and sine transforms
  yields a large class of fast general radix algorithms, many of which
  have not been found before.

Aside, when those transforms are considered in the filter bank framework, other fast versions exist, based on polyphase decomposition or lifting transforms.

Additional references: I do not know of a full book on  Algebraic Signal Processing Theory. The above link has a handful of references. For book-style, you have PhD theses:

A. Sandryhaila, 2010, Algebraic Signal Processing: Modeling and Subband Analysis
M. Püschel, 2008, DFT and FFT: An Algebraic View,  Computer Algebra Handbook, Foundations, Applications, Systems, Eds. J. Grabmeier, E. Kaltofen, V. Weispfenning

On filter banks:

L. Liu,  On Filter Bank and Transform Design with the Lifting Scheme
K. Soman et al., 2020, Insight into Wavelets: from Theory to Practice


How to find out whether guassian or salt& pepper noise present in an image, Is there any algorithm

Exploit the fact that salt&pepper noise has only 2 values and median filter does excellent job in getting rid of it.

Filter your image with median filter
Find pixels where median filter did most of its job
Plot histogram of those pixels
Compute fraction of 0s and 1s

threshold = 0.2; % adjust to perform best on your images
im = imread('../lena.bmp');
imd = double(im);
imd = imd./max(max(imd));
imdg = imnoise(imd,'gaussian');
imdsp = imnoise(imd,'salt & pepper');
subplot(2,2,1), imshow(imdg)
subplot(2,2,2), imshow(imdsp)
imdgf = medfilt2(imdg);
imdspf = medfilt2(imdsp);
g = imdg.-imdgf;
sp = imdsp.-imdspf;
gm = abs(g)>threshold;
spm = abs(sp)>threshold;
subplot(2,2,3), hist(imdg(gm))
subplot(2,2,4), hist(imdsp(spm))
waitforbuttonpress


How to create a band limited signal of sine waves that is nearly uniform distributed

This is a wierd thing you are trying to do. What i the purpose behind this signal generation? What particular application needs this?

I think of two ways, but maybe one is not correct.

1 - Generate Uniform White noise, and then filter it to your desired bandwith of interest. But that will surely affect the probability distribution of the amplitudes and make it more laplacian, like your graph shows.

2 - Generate Any signal you want with the desired narrow bandwith (Still not sure why you need this, but let's say you do). you could use method one for it. Then estimate the probabiliy density with any method, like periodogram like you did, or Kernel Density estimation with parzen windows, etc.
Then, design a non-linear transformation that you can apply sample to sample, in order to go from the probability density that you have, to a uniform. This non linear transofrmation is actually the "Inverse Probability distribution" that you have.

If you estimate your density in discrete intervals with periodogram, maybe you could look at the form of the graph and try to guess what functional form it has, then fit it with least squares to a continuous function, and then obtain the Probability Distribution from that, find the inverse, and apply that transformation to every sample. That should give you a signal which is uniformly distributed.

I am guessing now, but i think that your bandwith will not be preserved by this transformation. I have actually never studied the relationship between bandwith and probability density of the amplitudes. But try it and see what happens.

hope it helps.

Edit:
I think you are trying to do is Least Squares fitting. It is the most mathematically correct application of "averaging to get rid of interferences and unwanted noise", or of linearizing a non-linear systems. Take all your measurments, and fit them to a line. There are 5 types of Least Squares methods:

Least Squares with $L_{2}$-norm
Least Squares with $L_{1}$-norm
Least Squares with $L_{\infty}$-norm
Total Least Squares
Data Least Squares

Try some of them and see which one gets better result for your particular application.


one thing about a non-minimum phase system (with a rational transfer function), is that it can be thought of as the series concatenation (or cascade) of a minimum-phase system, having identical magnitude response as the given non-min-phase filter, with an all-pass filter.  the APF will have a poles that cancels specific zeros of the min-phase system that are inside the unit circle, and puts in new zeros ("reflected zeros") outside the unit circle.  so the new system always has more phase shift than the min-phase system.

so here's another definition of a minimum phase system:

An LTI system or filter

$$ H(f) = |H(f)| e^{j \arg\{H(f)\}} = |H(f)| e^{j \phi(f)} $$

is minimum phase if and only if the natural phase response, in radians, is the negative of the Hilbert transform of the natural logarithm of the magnitude response:

$$ \phi(f) \triangleq \arg\{ H(f) \} = -\mathscr{H}\big\{ \ln( |H(f)| ) \big\} $$

since

$$\begin{align}
 H(f) & = |H(f)| e^{j \phi(f)} \\
  & = e^{\ln(|H(f)|)} e^{j \phi(f)} \\
  & = e^{\ln(|H(f)|) + j \phi(f)} \\
  & = e^{\ln(H(f))} \\
\end{align}$$

this is relating the real and imaginary parts of the complex natural logarithm of the frequency response.

$$ \Im\{ \ln(H(f)) \} = -\mathscr{H}\big\{ \Re\{\ln(H(f)) \} \big\} $$

so a question to ask, why is the pole/zero definition stated in the question equivalent to this Hilbert transform definition??


Reconstructing time domain signal with Hanning window

You should be using overlap.  That ought to take care of the "oscillating" and make it unnecessary to undo the Hanning window.

The effect you are seeing comes about because the window values at the edges get very close to (and may even be) zero.  The window edges won't line up perfectly with the signal after you've done your processing in the frequency domain, so dividing by those very small values will cause the peaks you see.

In your place, I would first make sure that I could do the windowed FFT/IFFT and recover the original signal.  When you can do that without getting the "oscillation" you can put your processing back in.


The frequency domain conjugate multiplication (correlation) of the received signal with the reference signal followed by the power delay profile will provide you the overall signal to noise ratio. As explained in paper titled "SNR Estimation based on Sounding Reference Signal in LTE Uplink".

There are various noise reduction algorithm exist. Which will reduce the noise and will enhance the main expected signal component. There are adaptive filter which will converge to the expected value and depend on there performance, complexity is dependent. For example LMS, NLMS and RLS are the type of adaptive filters. Which will reduce the noise but the accuracy is directly proportional to the complexity.

By filtering we can reduce the noise, but we cant remove it totally for example AWGN noise.


The current ImageMagick can make them.  Given a little-endian TIFF,

convert le.tiff -endian msb -define tiff:endian=msb be.tiff

and vice versa:

convert be.tiff -endian lsb -define tiff:endian=lsb le.tiff

Run

identify -verbose be.tiff | grep ndian

to see what happened.


Note that

$$1+z^{-1}+\ldots + z^{-(N-1)}=\sum_{n=0}^{N-1}z^{-n}=\frac{1-z^{-N}}{1-z^{-1}}\tag{1}$$

where I've used the formula for a finite geometric series.

So both your results are identical and correct.

The ROC is $|z|>0$, which is the case for all causal sequences of finite length. Note that in the expression on the right-hand side of $(1)$ there is a pole-zero cancellation at $z=1$, so in fact there is no pole at $z=1$, hence the ROC $|z|>0$.


How to perform downconversion using a DSP

Yes, you can do the operation digitally; this is a very common operation in modern software-defined radio systems. The implementation is straightforward: you just implement the same steps that you would see in a hardware implementation. First, multiply by a sinusoid at the carrier frequency, then apply a lowpass filter to remove any images.


How to achieve an arbitrary phase shift?

You can't design a filter that creates a phase shift that's constant with frequency for  real valued input (if that's what you are trying to do).

A Hilbert transformer appears to be doing this. However, the problem is, you can't implement a perfect Hilbert transformer since it's non causal with an infinite length impulse response.

The tricky part is that the DC and Nyquist components of the spectrum have to be real valued and that the phase at these frequencies is constrained to $0$ or $\pi$. Any phase shifter would have to have a transition band from DC to the to target phase and from the target phase to Nyquist.

That's exactly what happens if you design a real-world Hilbert transformer: you trade off the size of the transition bands against the length & complexity of the filter.

The class of filters that only manipulate phase and not the amplitude are called allpass filters. One can easily show that

Allpass filters have a phase that's monotonously decreasing
For an allpass filter of oder $N$ the phase from DC to Nyquist decreases by $N \cdot \pi$
Allpass have zeros that are inverses of the poles,
The z-transform numerator polynomial is the reverse of denominator polynomial

Of course you can always design an "approximation" that's good enough for your requirement. You'd have to specify the target bandwidth, max phase error, max amplitude error and then deploy a suitable FIR or IIR design method

Even an approximate design of absolute phase shift is quite difficult, especially if you have latency or causality constraints. Constant relative phase shift is a lot easier. A popular technique is to run a signal  through two parallel allpass filters which are designed to have a constant phase difference over the frequency area of interest. The outputs have identical amplitude as the input but have approximately constant phase shift with respect to each other (although not with respect to the input)


It depends. Theoretically you can save a few multiplies with a radix-4 as radix-4 has a 1/4th the number of butterflies and 3 mpy + 8 adds per butterfly (if properly structured) and the radix 2 has 1 mpy + 2 adds per butterfly.

So in terms of multiplies it's a bit better, however there is higher complexity in terms of code structure, exception handling, coefficient management, register management, digit-reverse addressing, etc.

So it's only an advantage if the number of mpy is the limiting factor which for most hardware these days is not the case.


Matching LPC frequency response to raw sample data frequency response

Is using the coefficients directly as input to the FFT the right approach?

It depends on what are you going to do next. There are many formant extraction algorithms, you can find recent advanced one in COVAREP https://github.com/covarep/covarep/blob/master/formant/formant_CGDZP.m

  and perform FFT directly on the LP coefficients padded with zeroes.

It seems like a crazy step with not much meaning. Why are you doing this? LP coefficients create very rough spectrum approximation, the rest is in residual, but FFT from LPC, it is something strange.

Windowing does not change much, the difference will not be visible. You only can get a few percent more accurate results from it.


RC high pass filter differential equation

In an RC high-pass filter, the output is taken across the resistor instead of the capacitor. Otherwise, if the output was taken across the capacitor, the RC high-pass filter would start to work as an RC low-pass filter.
An example RC high-pass filter circuit is given below:

The basic idea of using an RC high-pass filter is to filter out the input signal’s low frequency components and let the high frequency components to pass through the succeeding stages of a network.
Actually, this filtering operation depends heavily on the capacitor’s properties for varying frequencies. The formula of a capacitor’s capacitive reactance which represents its resistive nature for varying signal frequencies is given below:
$$X_C{(j\omega)} = \frac{1}{j\omega{C}} \Omega$$
As it can be seen, the formula results in a complex value. According to the formula above; if the frequency $f$ of the signal which is passing through the capacitor decreases, so do the angular frequency $\omega$ because of the formula $\omega = 2\pi{f}$, the magnitude of the capacitive reactance increases. Therefore, basically, the signal which is passing through the capacitor is attenuated.
So, in the light of these knowledge, let’s analyse the given circuit.
If the Kirchhoff’s voltage law is applied on the circuit, we get;
$$-V(t) + V_C{(t)} + V_R{(t)} = 0$$
$$-V(t) + {\frac{1}{C}\int_{t_0}^{t} i(t) \, dt} + V_C{(t_0)} + V_O{(t)} = 0$$
If the assignments $V(t) = x(t)$ and $V_O{(t)} = y(t)$ are done in order to involve classical input and output notations as in the systems, the equation becomes;
$$-x(t) + {\frac{1}{RC}\int_{t_0}^{t} y(t) \, dt} + y(t) = 0$$
$${\frac{1}{RC}\int_{t_0}^{t} y(t) \, dt} + y(t) = x(t)$$
$${\int_{t_0}^{t} y(t) \, dt} + RCy(t) = RCx(t)$$
If the time derivatives of the both sides of the last equation is taken, the resultant equation becomes;
$$y(t) + RC\frac{dy(t)}{dt} = RC\frac{dx(t)}{dt}$$
According to this differential equation which describes the input-output relationship of the given RC high-pass filter, the frequency domain analyses can be done.
The same process can be followed when analysing an RC low-pass filter too.


How much can the FFT noise affect the timbre analysis of a harmonic sound?

An FFT is a finite length transform.  Therefore one does need to take into account any artifacts introduced by the window (implied rectangular, or optional Von Hann, et.al.).

There is some controversy, but you might consider sampling using a greater ADC bit depth than 16 bits at a rate higher than 44.1 ksps to feed you FFTs, and use one of the typical double precision FFT libraries to reduce any numerical noise.  (...and how you mic the violin may make a far greater difference than any of the above.)


extracting MPEG-7 low-level audio descriptors using essentia

Looking at the Essentia documentation, the MPEG-7 isn't listed directly as a function from where you can extract features.

The audio descriptors in the MPEG-7 standard are listed here. I think you will have to take use various functions to create your own feature vector. Hope this helps!


A classical result states that a (non-trivial)  time-limited signal does not have a  limited bandwidth (in terms of support). The converse happens as well: a bandwidth-limited signal cannot be of finite length.

And there is a similar rule-of-thumb principle:  when you try to reduce the support of a signal in one domain (time or frequency), you generally spread it in the other domain.


If you have the blurred and unblurred images, what I think you are asking is how you can recover the point spread function (PSF).

In the absence of noise, this is theoretically possible by considering the operation in the frequency domain. The blurred was introduced by a convolution of the unblurred image with a Gaussian kernel. In Fourier space, convolution is applied by multiplication of the image with the kernel, so if you have two of the operands (the source and the result), it should be possible to recover the blur kernel by division of the respective Fourier transforms.

This is effectively what some deconvolution operations attempt to do when the PSF is unknown or is a mixture of defocus and motion blur (e.g. blind deconvolution, Wiener filter, etc.).

So, in the frequency domain, the blur can be expressed as

G(u,v) = H(u,v)F(u,v) + N(u,v)

Where G is the Fourier transform of the blurred image, H is the blurring function, F is the source image and N is additive noise. Assuming the noise is negligible (unrealistic in most real world situations), the blur kernel can be recovered by calculating

H(u,v) = G(u,v) / F(u,v)

and then taking the inverse Fourier transform of H to obtain the PSF as an intensity image. From this you could measure the radius.


Does a PID need to be "online" to be tested? Can the output signal be viewed without it being used?

The PID can be tested in an open-loop setup, however any integral gain <> 0 will cause the PID output to run away, so you are likely to hit clamping (like control voltage maxima) with your test rig.

What really matters and needs to be tuned is the closed-loop behaviour. This includes your plants behaviour. Stability is prime goal, maximum overshoot and settle time (to set-point change or disturbance) can be tuned.

If you cannot tune in the live system, modeling and simulation is the way to go (look for Simulink Control Design). You would typically start with an empirical plant model, tune your PID for it, then run the actual device with this PID but a set-point well below Tmax. Instrumentation will give you data to improve your plant model. Repeat this development cycle a few times until you are confident to increase to the actual set-point.


In a wireless system, most of the noise is added by the receiver itself. The received signal is very weak (-100 dB or less is typical), and its power is comparable to the thermal noise of the receiver's amplifiers. That is the reason why noise is modeled as being added to the signal present in the antenna.

The noise still affects the estimated signal, of course; if the received vector is $r$, a naive ZF receiver would calculate

\begin{align} \hat{y} &= H^{-1}r \\ &= H^{-1}[Hy+n] \\ &= y + Hn. \end{align}

Since the noise $Hn$ is correlated, this receiver has terrible performance.


how sinc pulses can be used to distinguish different symbols

How can sinc pulses be used to distinguish the different numbers? If we are using 16-QAM, then we have 16 different symbols. Is it the amplitude of the sinc pulse that will change to distinguish the 16 symbols?

You send one sinc after the other, and each sinc is multiplied with one complex number (which is your 16QAM constellation point)

The bandwidth of a single sinc pulse is Rp/2, what is the bandwidth of the whole signal s(t) consisting of several pulses, each of width Rp/2?

The bandwidth of something is something that can be seen from the power spectral density, which in itself is an average.
Hence, the spectrum of your transmission, assuming uncorrelated data symbols, is exactly the spectrum of your pulse.

Your two questions are the most fundamental statements underlying all modern digital communications:

you can transport data by sending a symbol impulse train through your pulse shaping filter
the spectrum of the signal is defined only by pulse shape

It feels a bit like you learned about 16QAM before learning the basics of complex baseband and modulation – maybe go a chapter back in your book until you find a block diagram showing a data source, a symbol mapper, a pulse shaper, and a mixer, and then refresh your knowledge of complex baseband, digital comms, if you feel like you're missing basics.


Finding repeating but not periodically repeating patterns of an unknown template which you expect the algorithm to identify is a hard problem.  Nature typically doesn't shut off all the other patterns that may be present for our convienance either.

Signals have sources or generating mechanisms. It would help if you had some idea how your patterns are generated.

The concept of noise can be helpful because what isn't noise might very well be signal.

I would approach a problem as a time series that corresponds to an ODE and try to learn the ODE.

One way to do this would be to construct a phase space representation, that is something like $x(t)$ and $y(t)=dx(t)/dt$. If the noise is low using central differences can give a good approximation to the derivative of your time series. A phase space plot might identify orbits and limit cycles produced by your pattern.

These ideas are common in chaos theory and I would search Google accordingly

As an example https://en.m.wikipedia.org/wiki/Van_der_Pol_oscillator


How can I infer the cost function from Kruppa's simplified equations

According to the equation(30) in this paper The solution was to reform the Equations from the picture in my question here, and after that(Which the paper mentioned in different, complicated mathematical format) is to take the Frobenius norm to the subtract of newly formulated equations(30) on the paper.

$\rho_{1}/ \phi_{1} = \rho_{2}/ \phi_{2} = \rho_{2}/ \phi_{2}$ where $\rho$ is the numerator of each term of the equations, and $\phi$ is the denominator of each term, now the equation(30) in the paper proposed that: $\rho_{1}* \phi_{2} - \phi_{1}*  \rho_{2} = 0 $ and $\rho_{1}* \phi_{3} - \phi_{1}*  \rho_{3} = 0 $ now taking the Frobenius norm $\parallel.\parallel_{fro} $ for the two equations yields two values, Thus the cost function is the sum of this values.


Technique of sampling a high-frequency signal (Mhz) with a lower sampling rate (khz)?

You cannot capture a 5 MHz signal by sampling at 48 kHz. That doesn't mean there won't be anything in your sample, it means it will represent the signal incorrectly.

What happens if you use a too small sampling frequency is called aliasing: the original frequency shows up as another frequency altogether. The translation from the real to the sampled frequency is like a process of wrapping or folding the infinite frequency axis into your limited frequency interval $[0 ; \frac12 f_s]$, or $[-\frac12 f_s ; \frac12 f_s]$ for complex-valued signals. $\frac12 f_s$ is also called Nyquist frequency.

In your example, not just the main frequency but the total bandwidth of the signal of 800 kHz is much larger than the sampling frequency. Because of this, I wouldn't expect any clearly defined peak to be present. Since you seem to find peaks, my guess is that that is due to your very primitive approach to spectral estimation. If you use a better method, e.g. Welch's modified periodogram, which is implemented in Matlab as pwelch, I'd expect that you see a very flat (white) spectrum.

Addition based on nimrodm's comment: My initial statement has to be relativized if the bandwidth of the signal is smaller than $\frac12 f_s$. In that case, if there is additional knowledge about the true signal frequency, this knowledge can be used to resolve the ambiguity introduced by aliasing, and perfect reconstruction of the true signal is possible. See undersampling. However, it is not possible to determine the true frequency from the undersampled signal and nothing else.

A simulated example: Matlab code to generate a frequency-modulated signal with peak at 5 MHz and bandwidth 400 kHz, initially sampled at 24 MHz.

% signal properties
f0 = 5e6;               % frequency 5 MHz
sf = 0.4e6;             % bandwidth 400 kHz

% sample properties
fs = 24e6;              % sampling frequency 24 MHz
N = 10e6;               % number of samples

% generate random frequency modulated sinusoidal signal
fi = smooth(randn(N, 1), 11);
fi = fi / std(fi) * sf + f0;
x = cos(2 * pi * cumsum(fi) / fs);

Estimating the spectrum using pwelch

% estimate spectrum
[psd, f] = pwelch(x, 512, [], [], fs);
figure
plot(f, psd)
xlabel('f [Hz]')
ylabel('normalized power spectral density [1/Hz]')

gives this result:

If the signal is downsampled to 48 kHz

% downsample to 48 kHz (factor 500)
fs = fs / 500;
x = x(1 :500: end);

the estimated spectrum looks like this:

This is, as expected, basically a white (flat) spectrum, and the many small "peaks" are random fluctuations due to the fact that this spectrum was estimated based on only a 500th of the original 10 million samples. The estimate could be improved by repeating the procedure and averaging over the different results.

If the bandwidth is changed to 4 kHz

sf = 4e3;               % bandwidth 4 kHz

the results of spectal estimation on original and downsampled signal look like this:

The narrower peak fits nicely into the range 0 to 24 kHz of the downsampled spectral estimate, it just necessarily appears at the wrong place, in this case 8 kHz.

8 kHz = 5 MHz – 104 · 48 KHz is the smallest alias of 5 MHz with respect to a sampling rate of 48 kHz, the next largest is 48 kHz – 8 kHz = 40 kHz (mirroring around the Nyquist frequency of 24 kHz). All the larger aliases derive from these two by adding integer multiples of 48 kHz: 56, 88, 104, 136 kHz up to 4952, 4984, 5000, 5032 kHz and so on.

Since the spectral estimation cannot distinguish these aliases, the frequency axis of the spectrum plot becomes ambiguous: instead of representing frequency components of the signal from 0 to 24 kHz, it could as well represent frequencies from 48 to 24 kHz (backwards because of the mirroring!),

or from 48 to 72 kHz, and so on.

If we know the frequency content of the signal lies around 5 MHz, we can select the corresponding aliased frequency interval, 4992 to 5016 kHz,

and recover the correct location of the peak at 5 GHz.

This trick does not work, however, if the bandwidth of the signal is larger than $\frac12 f_s$, because the true frequency content of the signal has been aliased across several of these frequency intervals.


Digital filter not removing noise at specific frequency Matlab

As mentioned in my comment, the filter returned by iirnotch is useless. From your filter coefficients you can see that the filter is only marginally stable due to two poles on the unit circle at DC and at Nyquist. Furthermore, even though the filter has a notch, it also attenuates all other frequencies quite strongly (apart from DC and Nyquist). The reason for that behavior is the extremely large bandwidth in your specification.

The figure below shows the magnitude responses of the filter you designed (top) and of a notch filter with a bandwidth BW = w0/35 (bottom) (note that the extremely large values very close to DC and Nyquist due to the poles are not shown in the top figure):

In any case, the bottom figure is what a notch filter should look like. If you tried that filter and it didn't do what you expected it to do, then the reason might be that your estimation of the noise frequency is wrong. Could it be that you got it wrong by a factor of $2$ (i.e., it would be a 60Hz hum)? [Also, doesn't the file name ecg58... suggest a disturbance at $2\cdot 29=58$Hz?]

So there might be several problems in your approach, but one is definitely the design of the notch filter, and if I may guess I would say that the other is the estimation of the noise frequency.


Relation between order and stability in IIR filter

Stability depends on the pole locations, so one cannot generally say that increasing the filter order does or doesn't affect stability. It depends on how you increase the filter order. What are the additional coefficients of the denominator? (Because the denominator polynomial of the transfer function determines the pole locations).

Another question are the finite word length effects. If you quantize the filter coefficients, e.g., by using fixed point arithmetic (but of course also by using floating point), then higher filter orders tend to be more problematic than lower filter orders, because the chance that poles get very close to (or even outside) the unit circle through quantization becomes larger. But again, you cannot say anything in general unless you specify the filter's transfer function, the filter structure, and the type of quantization.


Classifing audio signals to detect fault

You might want to look at:

Detection of Abrupt Changes:
Theory and Application1
Michele Basseville `
IRISA/CNRS
Rennes, France
Igor V. Nikiforov
Institute of Control Sciences
Moscow, Russia

with the free pdf at:

ftp://ftp.irisa.fr/local/as/mb/k11.pdf

It's a little old, but was advanced when it came out.

There are a lot of Google results for this topic.


How to write lowpass filter for sampled signal in Python?

The frequency response for the filter designed using the butter
function is:

But there is no reason to limit the filter to a constant monotonic
filter design.  If you desire a higher attenuation in the stopband
and steeper transition band, other options exist.  For more information
on specifying a filter using iirdesing see this.  As shown by
the frequency response plots for the butter design the cutoff
frequency (-3dB point) is far from the goal.  This can be alleviated
by down-sampling before filtering (the design
functions will have a difficult time with such a narrow filter, 2% of
the bandwidth).  Lets look at filtering the original sample
rate with the cutoff specified.

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

from scipy.signal import fir_filter_design as ffd
from scipy.signal import filter_design as ifd

# setup some of the required parameters
Fs = 1e9           # sample-rate defined in the question, down-sampled

# remez (fir) design arguements
Fpass = 10e6       # passband edge
Fstop = 11.1e6     # stopband edge, transition band 100kHz
Wp = Fpass/(Fs)    # pass normalized frequency
Ws = Fstop/(Fs)    # stop normalized frequency

# iirdesign agruements
Wip = (Fpass)/(Fs/2)
Wis = (Fstop+1e6)/(Fs/2)
Rp = 1             # passband ripple
As = 42            # stopband attenuation

# Create a FIR filter, the remez function takes a list of
# "bands" and the amplitude for each band.
taps = 4096
br = ffd.remez(taps, [0, Wp, Ws, .5], [1,0], maxiter=10000)

# The iirdesign takes passband, stopband, passband ripple,
# and stop attenuation.
bc, ac = ifd.iirdesign(Wip, Wis, Rp, As, ftype='ellip')
bb, ab = ifd.iirdesign(Wip, Wis, Rp, As, ftype='cheby2')

As mentioned, because we are trying to filter such a small percent
of the bandwidth the filter will not have a sharp cutoff.  In this
case, lowpass filter, we can reduce the bandwidth to get a better
looking filter.  The python/scipy.signal resample function can
be used to reduce the bandwidth.

Note the resample function will perform filtering to prevent
aliasing.  Prefiltering can also be perfomed (to reduce aliasing)
and in this case we could simply resample by 100 and be done, but
the question asked about creating filters.  For this example
we will downsample by 25 and create a new filter

R = 25;            # how much to down sample by
Fsr = Fs/25.       # down-sampled sample rate
xs = signal.resample(x, len(x)/25.)

If we update the design parameters for the FIR filter the
new response is.

# Down sampled version, create new filter and plot spectrum
R = 25.             # how much to down sample by
Fsr = Fs/R          # down-sampled sample rate
Fstop = 11.1e6      # modified stopband
Wp = Fpass/(Fsr)    # pass normalized frequency
Ws = Fstop/(Fsr)    # stop normalized frequency
taps = 256
br = ffd.remez(taps, [0, Wp, Ws, .5], [1,0], maxiter=10000)

The filter operating on the downsampled data has a better
response.  Another benefit of using a FIR filter is that you
will have linear phase response.


How to compute an old style shore line drawing?

Define a level function curve (L) like this:

$$ L(x,y) = C_0 \cdot x + C_1 \cdot y + C_2 \cdot x^2 + C_3 \cdot xy + C_4 \cdot y^2 .... $$

It has no constanst term.  Now, find the best fit to the level curve $L(x,y)=1$.  This will yield the values for the coefficients.

Normalize the coefficients so the gradient is of a magnitude you like.

Find the level curves of $L(x,y)=t$ where $t$ is near one.

Finding the best fit is done with Linear Algebra.  Finding the level curves is an interative search algorithm.


How to remove spike from the data series?

I would put this in a comment but I don't have enough rep points, as just started on this.

I would recommend reading this information before designing with sigma-delta ADCs
http://www.numerix-dsp.com/appsnotes/APR8-sigma-delta.pdf

It may be that you need a basic analog filter before the ADC input to prevent interfrence noise aliasing back in to your signal (assuming the load on your cell is constant?).


do i have an image processor to compare an image with a number of images and choose the identical image??

This is a very broad question and to be honest there are a multitude of ways to go about this. I highly recommend that you read more into imaging processing literature in general. From this you could perhaps suggest an approach that we could base our answers on, and help hone in on a good method for you.

That being said, a common approach to your problem is the SIFT algorithm. This algorithm will give back a set of features for a given image. From this you could construct a way to recognise images. However, this might not be the simplest solution if you are just a beginner at this game.

Best of luck!


The complexity of such function run in Matlab

You should know what each operator (i.e., *) and each called function (i.e. product) does.  Then add up those operations.  For instance, I'm pretty sure that in Matlab, the way you're building up Aug_t means that Aug_t' * Aug_t generates a vector dot-product, so for each element in Aug_t there's a multiply-accumulate.


Does waveform-based ASR with lowpass filtering preprocessing still require the assumption of FFT?

Pre-emphasis and low-pass filtering can be done completely in the time domain (using IIR or FIR filtering, for instance), thus without an FFT (or with, by fast convolutions filtering methods).


Transmission frquency in MIMO-OFDM system

For a given subcarrier, QAM symbols are transmitted on several antennas (MIMO).
And an antenna can transmit several subcarriers (OFDM). One can imagine that QAM symbols are modulated on a 3-dimensional grid frequency-space-time.


FFT strong signal in middle frequency bin

If your system is powered by a switching power supply that is not synchronized to the PRF (or multiples of), then you may get reliable spurs over time as seen in your histogram.

In your case, it's been found that was indeed the problem! The 500 KHz frequency fell nicely in the middle of your 1 MHz histogram. Hopefully this will help hunt down similar issues in the future.


The book uses the convention that $Z$ represents a unit delay (this is explained on this page). The more common convention is to denote the unit delay by $z^{-1}$. Consequently, minimum-phase systems have all their zeros inside the unit circle of the $z$-plane, but outside the unit circle of the $Z$-plane.


How critical is the selection of the window function in STFTs?

Window functions have an inherent tradeoff between two of their frequency-domain properties:

Main lobe width: Any tapered window function will cause some "smearing" in the frequency domain. This is visualized by the width of the center lobe in the window function's frequency response. The wider the main lobe, the more difficult it is to resolve two tones that are close in frequency (if they are closer to one another than the main lobe width, they will tend to smear together). So ideally, you would like to have a window function that has a very narrow main lobe.
Maximum sidelobe height: Many window functions have frequency responses that consist of a single main lobe surrounded by repeated sidelobes that decay at some window-specific rate. The height of these sidelobes can make it difficult to resolve two tones that are separated in frequency, but differ greatly in amplitude. So ideally, you would like to have a window function that has very low sidelobes.

The problem: if you decrease the main lobe width of a window function, the sidelobes will grow, and vice versa. So, you need to strike an application-specific balance when choosing a window, based upon the distances in frequency and amplitude that you expect between your signals of interest. Given specific parameters of your system, it's possible to choose a window that (hopefully) meets your requirements.

As far as choosing the length of your window (which is equivalent to choosing the length of the DFT), you're best served with making your observation as long as possible within the constraints that your application might impose (e.g. latency requirements, how long the signals of interest can be considered stationary, computational resources, etc.). Your ability to resolve in frequency is directly proportional to the observation length (measured in time, not necessarily based on the FFT length, which can be zero-padded with no improvement in frequency resolution).


Can the instantaneous frequency be always derived from an analytic signal?

This is because the phase term is changing with time at the rate of $t^2$ for the frequency ramp, and since the instantanous frequency is the derivative of the phase your result is as expected.

$m(t) = cos(2\pi f(t) t)$

$\phi = 2\pi f(t) t$

$f(t) = Kt$ for a linear frequency ramp

thus phase versus time changes at the rate of $t^2$

$\phi = 2\pi Kt^2$

and the derivative of the phase (which is the instantaneous radian frequency) is:

$\frac{d\phi}{dt} = 4\pi Kt$

(Divide this by $2\pi$ to get instaneous frequency, but the source of the 2x discrepancy should now be very clear).

See the answer to this DSP Challenge that I posted that further explains this:

Simulation of a Frequency ramp


Evaluation of Autocorrelation and Power Spectral Density of white noise through a filter

White noise is modelled as a process having spectral density $G_n(f) = k$, a constant.

After the filter, the resulting process $y$ will have spectral density $G_y(f) = G_n(f) |H(f)|^2$. And the autocorrelation will be the inverse Fourier transform of the spectral density.

This is a general result for any stable filter and finite power, wide-sense stationary input process (in this case we are pushing the boundaries because white noise has infinite power as defined above).


I think I've managed to implement your system, but in R rather than matlab (sorry, I have no access to matlab and scilab stopped working for me on OSX).

The implementation gives the following estimates of the position and velocity, and seems to work even with the large velocity noise.

The top plot is the true (red), noisy (black), and estimated (blue) position. The bottom plot is the true (red), noisy (black), and estimated (blue) velocity.  I've zoomed in on the bottom plot to see the true and estimated values because the noise is so large.

To your questions:

  My question is: how can I find Q and R?

Well, from my code Q is

  Q <- matrix(c(0,0,0,0,0,0,0,0,1),3,3)

which would be

  Q = [0,0,0; 0,0,0 ; 0,0,1];

in matlab, I think.

And R is

  R <- matrix(c(0.001,0,0,20),2,2)

which is

  R = [0.001 0; 0 20];

in matlab.

I've generated my acceleration as a vector that is T long with a variance of 1; you seem to have generated just one value for it with a variance of 100.

As for the second part of your question:

  Does it depend of what I want to do?

Generally, the better the values the Kalman filter use match the "true" values, the better the Kalman filter estimates the state.  I say "true" instead of true because sometimes we can't know what the truth is, so we have to guess it.

The only leeway I see in what you've presented is what the value if Q is. It's not clear to me what variance you assumed the acceleration to have.  If it is a constant, then the variance will be zero; if it's a Gaussian random variable, then you'll need to set the Q matrix using the variance of that random variable.

Also, this all assumes that the process noise and the measurement noise are uncorrelated (or independent). If there is correlation, then you'll somehow need to take that into account.

R Code

  # 26818
  T <- 128;
  dt <- 0.01

  t <- seq(0,T-1)*dt

  xkm1km1 <- matrix(c(0, 0, 0),3,1)
  Pkm1km1 <- matrix(c(1000,0,0 ,0,1000,0, 0,0,1000),3,3)
  #H = [1 0 0
  #     0 1 0];
  H <- matrix(c(1,0,0,1,0,0),c(2,3))

  #F = [ 1 dt dt^2/2
  #      0  1 dt
  #      0  0 1];
  A <- t(matrix(c(1,dt,dt^2/2,0,1,dt,0,0,0),3,3))
  B <- matrix(c(0,0,1),c(3,1))

  true_acceleration <- rnorm(T,0,1)
  measurement <- t(matrix(rep(0,T*2),c(2,T)))
  x <- t(matrix(rep(0,T*3),c(3,T)))

  position_noise <- rnorm(T,0,0.001)
  velocity_noise <- rnorm(T,0,20)

  v <- t(array(c(position_noise, velocity_noise),c(T,2)))

  Q <- matrix(c(0,0,0,0,0,0,0,0,1),3,3)
  R <- matrix(c(0.001,0,0,20),2,2)

  # Generate the data using the signal model
  for (k in 2:T)
  {
    x[k,] = A %*% x[k-1,] + B * true_acceleration[k-1]
    measurement[k,] = H %*% x[k,] + v[,k]
  }

  library("MASS") # For pseudo inverse ginv()

  zhat <- matrix(rep(0,2*T),c(T,2))

  for (k in 1:T)
  {
    xkkm1 <- A %*% xkm1km1
    Pkkm1 <- A %*% Pkm1km1 %*% t(A) + Q
    K <- Pkkm1 %*% t(H) %*% ginv( H %*% Pkkm1 %*% t(H) + R)
    z <- matrix(c(measurement[k,]), 2, 1)
    xkm1km1 <- xkkm1 + K %*% (z - H %*% xkkm1)
    Pkm1km1 <- (matrix(c(1,0,0,0,1,0,0,0,1),3,3) - K %*% H) %*% Pkkm1
    zhat[k,] <- as.numeric(H %*% xkkm1)
  }

  par(mfrow=c(2,1))
  plot(t, measurement[,1], type="b", pch=19)
  lines(t, x[,1], col="red",  lwd=10)
  lines(t, zhat[,1], col="blue", lwd=5)

  plot(t, measurement[,2], type="b", pch=19, ylim=c(-1,1))
  lines(t, x[,2], col="red",  lwd=10)
  lines(t, zhat[,2], col="blue", lwd=5)

  errs <- c(sum((signal - zhat)^2), sum((signal - lpf)^2))
  print(errs)


Optimal FEC option for extremely low SNR situation

The one I am considering implementing is an LDPC as outer code and a Turbocode as inner code, with an interleaver in between.

I have not seen any implementation of LDPC and Turbo Codes concatenation yet. There is probably a reason no one tried (which I don't know). The first limitation I see is the fact that decoders need soft-decision inputs in order to be optimum or close to theory. Building an LDPC or Turbo decoder with soft outputs is an effort (COST IPs deliver hard bits on outputs).
Inner-Outer coders topologies were used for technological reasons. At the beginning of communication systems main coders/decoders that were possible to implement (no IC technology limitations) were Convolutional Codes and Reed-Solomon. The decoding gain of each was not enough, so they were combined.
Then arrived Turbo Codes (in the 90's - I was a child) allowing more decoding gain in a single decoder than CC+RS. After then, in the 2000's arrived DVB-S2 and LDPC's. The reason why LDPC is combined with BCH is that there is an error floor. BCH is there to correct those last single errors LDPC cannot correct.

How do these codes perform against burst errors?

Interleaving is really to investigate case by case and depends on your channel. If you find that a burst of errors which can be to much to be corrected in a frame could happen, then interleave several frames is a possible approach.

Would any other combination of codes reduce the necessary Eb/No in order to achieve relatively low (<1e-04) BER?

I think you should be able to pick a single code. The choice is up to your link budget. Finding the code for a given architecture is another story important to consider!
There is a good document that gives an overview of codes used in space communications. Look for CCSDS 130.1-G-3. Figure 3-3 shows Shannon limits, Figure 3-5 some codes comparisons and Figure 7-10 shows Turbo Codes performances.
Also this website has an interactive FEC comparator: https://aff3ct.github.io/comparator.html
Whatever the code you choose to go to the lowest SNR or EbN0, you have to keep in mind how your demodulator is going to behave to lock the low signal. A 0 dB EbN0 LDPC 1/2 BPSK channel can be so "noisy" that some demodulators are not able to lock on it


Linearity and time-shifting of $\mathcal{F}\{0.8^n\cos(0.1πn)u[n]\}$

$\displaystyle \hbox{ DTFT}(\alpha x_1 + \beta x_2) = \alpha\cdot \hbox{ DTFT}(x_1) + \beta \cdot\hbox{ DTFT}(x_2)$

Linearity property cannot applied here.


Using Matched Filter To Predict the Existence of An Event

Yes. However cross correlating is matching so you got a mix-up in terminology. You can detect the presence of the signal this way.
Matched filter will always maximize the output SNR so it will not increase your noise. Given Gaussian noise it is also a Most Powerful test (I.e optimal).


Large spike at the center frequency when using Ettus X310

The problem that I have is that I always have a big spike (10-15 dB)
  directly on the center frequency (no matter what frequency I set). I
  am relatively new to all this so I would appreciate any pointers on
  how to get rid of the spike.

That spike is probably nothing surprising – just the LO leakage/DC offset, a very common artifact in direct conversion receivers!

Now, getting rid of it with a USRP is relatively simple.

You have to understand that the USRPs have a two-step tuning process:

There's a physical oscillator, the LO, which gets synthesized by the daughterboard, and is used to mix down the target signal to IQ baseband signals. Typically, synthesizers can only generate a discrete set of LO frequencies $f_\text{RF}$.
The ADC then samples these signals at the master clock rate, and the FPGA digitally shifts that sample stream in the master clock rate (i.e. 200 MHz for the X310 by default) wide baseband by essentially multiplying with $e^{j2\pi f_\text{offset} n}$, to give you your baseband signal from the desired $f_\text{target}$:

By default, $f_\text{RF}$ is set as closely as possible to $f_\text{target}$, minimizing $f_\text{offset}$, and hence, putting the LO into your observed bandwidth $b_\text{sample}$.

You can, however, manually set a "desired" $f_\text{offset}$; as soon as that is more than $\frac 12 b_\text{sample}$, the LO is out of your observer band!

Instead of using 101e6 in the Center Frequency input field, just use

uhd.tune_request(f_target, f_offset)

e.g.

uhd.tune_request(101e6, 4e6)

Let me comment on a few other issues directly, too:

  I am using the Ettus X310 SDR and GNU radio.

The X310 is the platform, indeed, but the actual analog signal processing (generating the LOs for mixing, mixing, amplification, baseband filtering) happens on the daughterboard. In such cases, it's usually important to specify which daughterboard you're using

  On the GNU radio schema I
  have a standard UHD USRP source block connected to a GUI FFT sink -
  basically I'm just displaying the spectrum of the raw signal from the
  device.

Indeed, that should work. Notice that the WX Gui blocks are kind of becoming deprecated; if you're just starting with GNU Radio and don't need to continue an existing application based on WX, I'd recommend changing the "build options" in the option block to "Qt", and replacing the "WX GUI FFT sink" with a "Qt GUI frequency sink". Same thing, more functionality, less CPU usage, usually.

  Parameters for the USRP source block:

  sample rate = 1.024 MS/s,

You probably were using a master clock rate of 184.32 MHz, right?

For future readers, reproducing this (this specific rate really comes up more often then you'd think!):

The X310 doesn't support that sampling rate in the default configuration! The console window should contain a clear warning that UHD decided to use a different frequency, and which frequency that was – in your case, probably 1.02048… MS/s.

UHD Warning:
    The hardware does not support the requested RX sample rate:
    Target sample rate: 1.024000 MSps
    Actual sample rate: 1.020408 MSps

That's because the X310 runs its ADC/DAC at what is called the master clock rate, which can take the values 200 MHz (default), 184.32 MHz (typical for some cellular network standards), 120 MHz (typical for LabView). The sample rate you deliver to the PC is a sample stream sampled at that rate, decimated in the FPGA to the rate you're requesting – and that decimation is done with "normal" decimating FIR filters, which can only give you an integer fraction of the input rate.

So,

$$1.020408… \frac{\text{MS}}{s} = \frac{200\frac{\text{ MS}}{s}}{196}$$

is possible,

$$1.024 \frac{\text{MS}}{s} = \frac{200\frac{\text{ MS}}{s}}{195.3125}$$

is not, when the MCR is set to it's default of 200 MHz.

With an MCR of $184.32\text{ MHz}$, which is $1.024\text{ MHz}*180$, this is indeed possible!

  center freq = 101 MHz,

Which probably means you're using the WBX or UBX daughterboard,

  gain = 42 dB,

which is not a gain setting that UBX or WBX support, their gain ranges are 0 – 31.5 dB and 0 – 31 dB, respectively.

  bandwidth = 130 MHz

That's the analog bandwidth specification; it applies to daughterboards with adjustable baseband low pass filters. However, there's no daughterboard in existence that has an adjustable baseband filter and a maximum bandwidth > 8 MHz.

In fact, this setting does nothing on UBX or WBX.


How to generate a noise like Pioneers white noise?

The "cool" stuff that you hear has almost nothing to do with the noise itself, but it's all in the post-processing: They use aggressive time variant low-pass and band-pass filtering, manual gain modulation and reverb/echo effects.
Any old pink noise will do just fine for this purpose.
"White noise" is an application specific misnomer. Almost no one uses actual white noise in audio, since it sounds really bad. Mostly it's pink noise or sometimes even brown noise.


Beat detection and calculation of start time

Interesting. I would be very tempted to use correlation and given the OP wants to specifically find the occurrence of a first pulse implies that the correlation interval is restrained to be one pulse cycle wide. (Unless we know further information, and specifically how many pulses occur after a start position- if this is known then further processing can be done on the entire "burst"). The first question to determine the applicability of an optimized correlation approach is how correlated each of the pulses are to each other, within a given sequence and from sequence to sequence (repeatability over time). This can be evaluated by cross-correlating various captures to a representative pulse and looking at the statistics on the peak magnitude of the cross correlation function, and how much margin this allows relative to the "noise" correlation everywhere else.  If that were to suggest correlation could be used (by providing a sufficient probability of detection vs probability of false alarm, which such metrics would be clear from this evaluation), I would then determine the mean reference pulse and use that to detect a pulse present or not present condition using a correlation function. This is basically and optimally a matched filter; optimally if the background noise is white. If the noise is not white, then targeted noise rejection filters could also provide significant advantages - this can be determined by capturing the background with no pulses present and evaluating the time and frequency statistics of those captures to see if there are particular details that can be targeted and eliminated (for example 60 Hz spurs from power lines, etc...).
Where a correlation / matched filter approach wouldn't work well is if the content of each pulse in time and frequency is mostly random from pulse to pulse, in which case an optimized receiver would instead consist of a power detector with a pre-selection filter around the dominant frequency content of the pulse prior to detection.  The first cross-correlation testing described above would reveal the optimum strategy based on those statistics. In either case, detecting the first (or any) occurrence limits the observation time for decision to be one pulse interval, which then limits the available SNR to make any such decision for any approach (unless other intelligence is known that make a first pulse unique from any other).


Noise floor estimation is usually done after applying an FFT to windowed data segments. By tracking the noise floor in each frequency band, the frequency dependence of the noise is taken into account. If the noise is non-stationary, its time dependence can be tracked by regularly updating the noise floor estimate in each frequency band. This usually requires the noise to vary more slowly than the desired signal. Have a look at this paper to get an idea of how it can work, and also look at the references therein. Search the web for "noise floor estimation" and "noise floor tracking" to find many more papers on the subject.


Ripples in a spectrum are an indication of echoes in time.  Specifically consider how the Fourier Transform of a sinewave is two impulses in frequency. Similarly, given the reciprocity of the Fourier Transform, a sinewave in frequency (in other words, ripples in the spectrum) is two impulses in the time domain. The Fourier Transform of the frequency response is the impulse response. A system that has an impulse response consisting of two impulses at different locations in time will output the input combined with a delayed copy of the input (echo).
In practice this can be caused by multipath delays for a signal propagating between a transmitter and receiver, reflections from poor matching on a high frequency (microwave) circuit, or other equivalent effects that result in multiple copies of the signal.


How to create a sine wave generator that can smoothly transition between frequencies

One approach that I have used in the past is to maintain a phase accumulator which is used as an index into a waveform lookup table. A phase delta value is added to the accumulator at each sample interval:

phase_index += phase_delta

To change frequency you change the phase delta that is added to the phase accumulator at each sample, e.g.

phase_delta = N * f / Fs

where:

phase_delta is the number of LUT samples to increment
freq is the desired output frequency
Fs is the sample rate

This guarantees that the output waveform is continuous even if you change phase_delta dynamically, e.g. for frequency changes, FM, etc.

For smoother changes in frequency (portamento) you can ramp the phase_delta value between its old value and new value over a suitable number of samples intervals rather than just changing it instantaneously.

Note that phase_index and phase_delta both have an integer and a fractional component, i.e. they need to be floating point or fixed point. The integer part of phase_index (modulo table size) is used as an index into the waveform LUT, and the fractional part may optionally be used for interpolation between adjacent LUT values for higher quality output and/or smaller LUT size.


The first equation, $ \frac{1}{N_s}\sum |x[n]|^2 $, is correct. Applying the RRC filter kind of smears the energy of the symbol over several samples so that when you add them all back up again you get $E_s$, and of course divide by $N_s$ for multiple symbols. The RRC filter should be normalized as someone in the comments noted.
The second equation I can't say much as it is not clear how you compute $P_{avg}$. Secondly, dimensional analysis is a fancy way of saying writing out the units and seeing what cancels out. For example, in the second equation:
$$\text{sps }T_s\rightarrow \frac{\text{samples}}{\text{symbol}}\frac{\text{seconds}}{\text{symbol}}=\frac{\text{(samples)(seconds)}}{(\text{symbol})^2}$$
I've found this method useful as a "sanity check" for myself many times. The above units should suggest to you that something is amiss (also why do you need two equations for the same quantity?).
Edit
The question was updated to include definition of $P_{avg}$. The OP is computing average power per sample and trying to convert it to average power per symbol. Again, write out the units it becomes clear:
$$\frac{\text{energy}}{\text{sample}}\frac{\text{samples}}{\text{symbol}}=\frac{\text{energy}}{\text{symbol}} $$
sps = 4;
Ns = 100;
Es = 10;
rrc = rcosdesign(0.25, 6, sps, 'sqrt');
s = sqrt(Es)*qammod(randi([0, 15], Ns, 1), 16, 'UnitAveragePower', true);
x = upfirdn(s, rrc, sps);

equation1 = sum(abs(x).^2)/Ns;
equation2 = sum(abs(x).^2)/length(x)*sps;


FFT Pitch Detection methods: Autocorrelation or other?

FFTs alone are lousy at pitch estimation for many really common sound sources, including male voices, and low piano or guitar notes.  Thus, the comments that an FFT is all you need are false and misleading.  (Although FFTs can be used as components of a more accurate composite pitch detection algorithm.)

When using autocorrelation, after you remove or discount subharmonic lag peaks, the pitch frequency may be around or close to the reciprocal of the time lag corresponding to the autocorrelation peak.

Pitch often refers to the human perception of such, and thus an estimate sometimes needs to include psycho-acoustics effects and illusions.  There are multiple methods for, and lots of book chapters and research papers on this problem of pitch detection.  How well each algorithm works might depend on your particular sound source and requirements.

Pitch detection/estimation methods include lag estimators, such as autocorrelation, weighted autocorrelation, AMDF and ASDF; and frequency domain analysis methods based on initial FFTs, such as cepstral/cepstrum methods, and the harmonic product spectrum algorithm.  Composite methods, and ones that use some decision analysis are less easy, but can be more accurate or robust for some sound sources.  Look for YIN, RAPT and YAAPT.  This list is by no means exhaustive.


Sound absorption coefficients for frequencies higher than 8kHz and lower than 80Hz?

The main reason why the tables don't have it, is that it's hard to measure. The measurement technique in ISO 354:2003 relies on measuring the difference in reverberation times in a reverberation rooms with and with/out a material sample.  At higher frequencies, the reverb time is dominated by air absorption and and the sound field becomes less and less diffuse which violates the basic assumption behind the measurement. It's not uncommon to measure absorption coefficients larger than 1. That's not good !

I wouldn't mingle air absorption and wall absorption. They are different effects and should be modelled separately. They are basically additive. At 20 kHz air absorption is around 50-ish dB/100m. The early reflections in a room have a travel distance of only a few meters so there is some attenuation, but there is still plenty of energy left.

Keep in mind that wall absorption measurement is NOT done at 8 kHz but it's averaged over the 8 kHz octave, so the data covers frequencies up to 11.5 kHz. If you need to go hire, you can often simply extrapolate to 16 kHz by looking at the difference between the 4 kHz and 8 kHz value.

Easiest is to extrapolate in dB for the reflection coefficient. Example: 4 kHz = 0.6, 8 kHz = 0.8. Reflection coefficients are 0.4 and 0.2 or -4dB and -7dB respectively. So it's 3 dB drop of reflected energy per octave. Extrapolation yields -10 dB for 16 kHz, which is a reflection coefficient of 0.1 or an absorption of 0.9.

Something like this
$$\alpha_{16} = 1-10^{0.1 \cdot (2 \cdot 10 \cdot \log_{10}(1-\alpha_8)-10 \cdot \log_{10}(1-\alpha_4))}$$


Is there any way to do 3D (Voxel Data) feature detection?

I never personally tried it but using image registration algorithm on Voxel Data you may have some results by trying registration on sub-cluster of voxel.
I would be interested in any feedback about this idea。


Calculation of actual analog input from bipolar ADC's output

Your understanding is right. First subtract the bias, and then apply scaling to get analog actuality.


This can't be done as a convolution because some of the operations you mentioned are not linear. Additionally, convolution by definition (multiplication by a toplitz matrix) shouldn't ever change the input and output size.

For bicubic:

http://en.wikipedia.org/wiki/Bicubic_interpolation#Bicubic_convolution_algorithm

Notice the high order terms.


In order for H(z) to be a linear phase filter, it must zeros both on the inside of the unit circle and at the complementary locations (1/z) which are outside the unit circle. Therefore a linear phase circle has no stable causal inverse (since this would necessitate poles outside of the unit circle.)

Note that a linear phase filter can be decomposed into the cascade of a minimum phase filter (all zeros inside the unit circle) and a maximum phase filter with the same magnitude response (all zeros outside the unit circle). For this reason, only causal minimum phase FIR filters have a stable causal inverse.

Further interesting detail, the reverse filter (all the coefficients in reverse order) of a minimum phase filter IS the maximum phase filter. This is all now intuitive: the coefficients for a minimum phase filter will be dominant toward the start of the filter (resulting in the minimum delay as the signal will emerge from the filter sooner). Reversing this results in the same magnitude response but the filter will have the maximum delay as the signal will emerge from the filter later. When you cascade two filters, you convolve their coefficients, so we can also see how the cascade of a linear phase filter with a maximum phase filter results in symmetric coefficients. And it is well understood that any FIR filter with symmetric coefficients is a linear phase filter!

Example [2 1 1] is a minimum phase filter, and [1 1 2] is the reverse maximum phase filter. The cascade of the two results in the convolution of [2 1 1] and [1 1 2] which is [2 3 6 3 2] which is a linear phase filter.

If the filter has real coefficients, the all zeros must also have a complex conjugate. (But whether the filter has real or complex coefficients does not change the answer--- both real and complex linear phase filters will not have a stable causal inverse, nor will they have a stable anti-causal inverse).

Regarding causality and the ROC, for causal FIR filter the ROC is $|z|>0$ (all the poles are at $z = 0$), and for anti-causal filters it is $|z| < \infty$ (and all the poles are at $z = \infty$). In either case the requirement for zeros to be at reciprocal locations holds in order for the result to have linear phase.


How does sampling and subsequent resampling affect aliased harmonics?

The alias behaves exactly the same, in every sense, as if it were an actual sinusoidal component at the alias frequency.

In your example, the 7 kHz harmonic sampled at 3 kHz is actually a 1 kHz sinusoid, and decimating it is exactly the same as decimating a 1 kHz signal. In this case, you end up with a DC value.


pixel data of color-images contain three (bytes of) intensity values denoted as R,G and B. Most common they are 8-bits but other formats are possible too.

Each intensity value at a pixel point represents the amount of light captured by a corresponding sensor at that point which is most sensitive to the given range of wavelengths (color); for example a RED value means that the associated RED sensor at that pixel captured that much of light in the wavelengths about $550$ nm to $750$ nm.

In a practical setting, the distinct wavelength captures can be realized by three different optical filters that separate the light spectrum into three regions associated with Red, Green and Blue. Then the incoming light (after being filtered) will be captured by three same type of sensors which are sensitive to the whole optical range (+ more possibly)


When is cubic spline interpolation better than an interpolating polynomial?

The phenomenon being discussed is Runge's phenomenon.

The maximum absolute value of the $n$th derivative of $\sin(\omega t)$ is $\omega^n$. For Runge's function $\frac{1}{25t^2+1}$ the maximum absolute value of the $n$th (even) derivative is $5^nn!,$ where $n!$ denotes factorial. This is much faster growth. Only if the derivatives grow too fast by increasing $n$, then it is possible that interpolation error diverges as the interpolation order is increased. Exponential in $n$ is not yet too fast. Have a look at: James F. Epperson, On the Runge example, The American Mathematical Monthly, vol. 94, 1987, pp. 329-341.

If a function has only continuous derivatives, then the competing approach, piece-wise polynomial spline interpolation always converges if a small fixed number of its early derivatives are bounded over the interval of interest, see Wikipedia article on linear interpolation as an example.

If both methods converge, then (non-piecewise) polynomial interpolation has the advantage of a higher polynomial degree if many samples are used, and may provide a better approximation, as you saw in your sine example. You may also be interested in L.N Trefethen, Two results on polynomial interpolation in equally spaced points, Journal of Approximation Theory
Volume 65, Issue 3, June 1991, Pages 247-260. Quote:

  [...] in band-limited interpolation of complex exponential functions
  $e^{i\alpha x}\, (\alpha \in \mathbb{R}),$ the error decreases to $0$ as $n \to \infty$ if and only if $\alpha$ is
  small enough to provide at least six points per wavelength.

You have 6.5 samples per wavelength.


For my DIY projects I have used ps3 eyetoy camera which can capture 120 fps. Only downside was difficulty in removing IR block filter with a filter blocking the visible range.
There used to be open drivers for this camera.
Hope it helps.


How to remove historical values from exponential filter

Yes. Your filter $f(t)$ is an LTI system, so you can use that to your advantage. Let:

$$
C(t) = B(t) - A(t)
$$

Then, if $y_A(t)$, $y_B(t)$, and $y_C(t)$ are the system's responses to $A(t)$, $B(t)$, and $C(t)$ as inputs, then you can straightforwardly do what you want:

$$
A(t) = B(t) - C(t) \rightarrow y_A(t) = y_B(t) - y_C(t)
$$

So, pass $B(t)$ through the filter, then pass $C(t) = [1, 2, 3, 0, 0, 0]$ through the filter and find the difference between the two filter outputs. That's exactly what you would get if you passed $A(t)$ through the filter on its own.


Quadrature demodulation, homodyne detection, lock-in detection - what's the difference?

They are essentially the same but may have key differences varying across implementation. The process is identical in terms of multiplying a signal by a sine wave at the center frequency of the signal (the carrier) which moves it to DC, and then low pass filtering that signal to maximize the SNR for the detection (presence of the signal at the carrier).

For the differences that may exist the homodyne and lock-in amplifier may or may not be done with quadrature demodulation. Note the following from the relationship of multiplying two sinusoids:

$$\cos(\alpha)\cos(\beta) = \tfrac{1}{2}\cos(\alpha+\beta) + \tfrac{1}{2}\cos(\alpha-\beta)$$

If we multiply two signals of the same frequency in phase and low pass filter to remove the higher frequency term (as denoted by the LPF[]) we get:

$$\mathrm{LPF}\big\{\cos(\omega t)\cos(\omega t + 0°)\big\} = \tfrac{1}{2}\cos(0) =\tfrac{1}{2}$$

but if in quadrature we would get:

$$\mathrm{LPF}\big\{\cos(\omega t)\cos(\omega t + 90°)\big\} = \tfrac{1}{2}\cos(90°) = 0$$

So a key question is what are you trying to detect? When trying to detect phase, the quadrature demodulator would actually be best as it would have the highest sensitivity to phase changes (the slope of the output voltage vs phase shift of the two signals being multiplied is highest when the signals are in quadrature). Similarly it will detect with maximum amplitude the presence of a signal that is also in quadrature (as this would be the first case essentially if both signals are at 90°). For arbitrary detection of any angle you would use an IQ demodulator that would multiply the input signal by cosine and sine with I and Q outputs each low pass filtered to show the detection for any phase of the input.


Applying 2D Image Convolution in Frequency Domain with Replicate Border Conditions in MATLAB

The trick is to properly compensate for the fact that Frequency Domain multiplication applies a convolution with the circular boundary conditions in the spatial domain.
You may use the following code:
clear('all');
close('all');

gaussianKernelStd       = 2;
gaussianKernelRadius    = ceil(5 * gaussianKernelStd);

mI = im2double(imread('cameraman.tif'));
mI = mI(:, :, 1);

numRows = size(mI, 1);
numCols = size(mI, 2);

vX = [-gaussianKernelRadius:gaussianKernelRadius].';
vK = exp(-(vX .* vX) ./ (2 * gaussianKernelStd * gaussianKernelStd));
mK = vK * vK.';
mK = mK ./ sum(mK(:)); %<! The Gaussian Kernel

mIPad = padarray(mI, [gaussianKernelRadius, gaussianKernelRadius], "replicate", 'both'); %<! For the replicate array

mKC = CircularExtension2D(mK, size(mIPad, 1), size(mIPad, 2)); %<! Circular extension for the 2D Kernel

startIdx = gaussianKernelRadius + 1;

mIFiltered      = ifft2(fft2(mIPad) .* fft2(mKC), 'symmetric');
mIFiltered      = mIFiltered(startIdx:(startIdx + numRows - 1), startIdx:(startIdx + numCols - 1)); %<! Removing the padding
mIFilteredRef   = imfilter(mI, mK, 'replicate', 'same', 'conv'); %<! Reference

figure(); imshow(mIFiltered);
figure(); imshow(mIFilteredRef);

max(abs(mIFilteredRef(:) - mIFiltered(:))) %<! Should be very very low

The function CircularExtension2D() is given In my StackExchange Signal Processing Q38542 GitHub Repository.  It was taken from my answer to Applying Image Filtering (Circular Convolution) in Frequency Domain.
The steps the code implements are as following:

Pad the image in order to have Replicate boundary condition convolution.
Convert the spatial domain kernel into a form which matches the image in frequency domain. We assume top left of the image is (0, 0) in spatial domain. So we need the (0, 0) of the kernel to also be in the top left corner.
Apply circular convolution using frequency domain.

As you can see, the result is prefect.
In my answer to How Much Zero Padding Do We Need to Perform Filtering in the Fourier Domain? I implemented a MATLAB Function, ImageFilteringFrequencyDomain(), to apply Frequency Domain Convolution with the border conditions supported in imfilter().
Remark
In your case, the kernel is 5x5 which is very small.
For small kernels and a single image it is better to apply the convolution in the spatial domain. It will be much faster.


Interpretation of FFT result when frequency is an exact multiple of frequency bin

Short term FFTs apply a windowing function at moving intervals, do an FFT, and average.  In this case, the interval with which the windowing function is applied is phase-locked to the signal, so what tends to be a process averaging out over the integration time no longer does.  You get the FFT of this intermediate product without averaging and it depends on just where you locked phase what comes out.

However, you did not state just what kind of signal you are using.  For a sine wave, I still would expect the result to be pretty close to a single peak.  It may be that this peak is absent due to fixed point overflow since the analyzer might not be prepared for the full concentration of frequency.  If that is the case, just reducing the amplitude might change the outcome to something more expected.


3D point cloud to 2.5D graycale image (metric) [Halcon?]

Your 3 channels already give you the 3D coordinates, and from that 3 images, you already know which pixels locations the 3D data are projected. So you don't need to do anything. If you like to get a depth image, just set:
$$
depth(i,j) = z(i,j)
$$
In Halcon, you just select the Z-channel using access_channel.
If you desire a distance image, then:
$$
distance(i,j) = \lVert (x,y,z)^T \rVert
$$
The norm defines the distance metric you would like to use. Traditionally, I would say Euclidean $\|\cdot\|=\sqrt{x^2+y^2+z^2}$. You could do this in any language by using operations like addition, multiplication and square roots.


Detecting a Signal Inside a Time Series Using MATLAB's $ \tt xcorr() $

The function xcorr calculates the correlation of 2 signals.
The correlation is known to be a good (The MLE) for delay estimation under Gaussian Noise.

Yet, as can be seen in your data you're not using it in the cases it meant to be used.
If we assume you have a model of a known signal with Additive White Gaussian Noise (AWGN or any other Additive White Noise), than the SNR is so terrible you can not expect it perform well.

If your signal model is additive deterministic signal (The signal [0 9 8 0 0 0 7 3 1 0] in that example) then you should look for a different approach (Though xcorr can be part of it).

Any tool you use, you must be aware for the model it was designed for.


If you would like to do this faster, what you could try is the neighborhood operations. Take for example the sample code below:

I = imread('tire.tif');
f = @(x) max(x(:));
I2 = nlfilter(I,[3 3],f);
imshow(I);
figure, imshow(I2);

This applies a maximum filtering. You could consider changing your neighborhood function f to your custom algorithm (inside your loop) and this should give you a better speed.

Of course there are a lot of ways of getting speed. Use of integral images, moving averages and precomputations might save you from the burden of heavy block processing. But this would be way a lot to write in here. For more details check out these papers:

http://www.csse.uwa.edu.au/~shafait/papers/Shafait-efficient-binarization-SPIE08.pdf
http://arxiv.org/ftp/arxiv/papers/1201/1201.5227.pdf

Cheers,


Real and complex low pass equivalent of band-pass signal

It's not true that any real-valued band pass signal can be written as shown in your first equation:

$$s(t)=A(t)\cos(\omega_c t)\tag{1}$$

The most general representation of a real-valued band pass signal is given by

$$s(t)=\Re\{a(t)e^{j\omega_c t}\}\tag{2}$$

where $a(t)$ is a complex-valued function, called the complex envelop of $s(t)$. It is the complex baseband (low pass) representation of $s(t)$. With $a(t)=a_I(t)+ja_Q(t)$, (2) can be written as

$$s(t)=a_I(t)\cos(\omega_c t)-a_Q(t)\sin(\omega_c t)\tag{3}$$

where $a_I(t)$ and $a_Q(t)$ are called the in-phase and the quadrature components, respectively. Both $a_I(t)$ and $a_Q(t)$ are real-valued baseband (low pass) signals.

A third, equivalent, representation of $s(t)$ is obtained by writing the complex envelop as

$$a(t)=b(t)e^{j\phi(t)}\tag{4}$$

where $b(t)$ is the (real-valued) envelop, and $\phi(t)$ is the phase. Using (4) and (2), the band pass signal $s(t)$ can be written as

$$s(t)=b(t)\cos(\omega_ct+\phi(t))\tag{5}$$

Equations (2), (3) and (5) are three equivalent representations of a general real-valued band pass signal. If you compare representation (5) to Equation (1), you see that Equation (1) is a special case for $\phi(t)=2\pi k$, $k=0,1,\ldots$, i.e. for a constant phase which is zero (or a multiple of $2\pi$). From Equation (3), this is equivalent with the condition that the quadrature component $a_Q(t)$ is zero.

Note that the complex-valued band pass signal

$$\tilde{s}(t)=a(t)e^{j\omega_c t}\tag{6}$$

is the analytic signal mentioned in a comment. It has no negative frequency components and it can be obtained from $s(t)$ by computing

$$\tilde{s}(t)=s(t)+j\mathcal{H}\{s(t)\}\tag{7}$$

where $\mathcal{H}\{\cdot\}$ denotes the Hilbert transform.


Determining Type and Bandwidth of a filter

As others have mentioned, performing a 2D FFT on the kernel will give you the frequency response of the filter. However, it's worth mentioning that 2D filters can be analyzed using the Z-transform, which may or may not provide deeper insight, depending on the filter (and what you want to know).

For example, given the kernel you specified, the corresponding difference equation would be

$$y(n_1,n_2) = x(n_1+1,n_2) + x(n_1,n_2 + 1) + x(n_1-1,n_2) + x(n_1,n_2-1) - 4x(n_1,n_2).$$

Its Z-transform is

$$ Y(z_1,z_2) = z_1X(z_1,z_2) + z_2X(z_1,z_2) + z^{-1}_1X(z_1,z_2) + z_2^{-1}X(z_1,z_2) - 4X(z_1,z_2),$$

which, after rearranging, yields the following transfer function for the filter:

$$ H(z_1,z_2) = \frac{Y(z_1,z_2)}{X(z_1,z_2)} = z_1 + z_2 + z_1^{-1} + z_2^{-1} - 4.$$

To determine the magnitude response, just plug in a pair of complex exponentials and simplify, as follows:

$$ \begin{aligned}
H(e^{iw_1},e^{iw_2}) & = e^{iw_1} + e^{iw_2} + e^{-iw_1} + e^{-iw_2} - 4  \\
                     & = (e^{iw_1} + e^{-iw_1}) + (e^{iw_2} + e^{-iw_2}) - 4 \\
                     & = 2\cos{w_1} + 2\cos{w_2} - 4 \\
|H(e^{iw_1},e^{iw_2})| & = 2\sqrt{ (\cos{w_1} + \cos{w_2} - 2)^2  } .
\end{aligned}
$$

Evaluating the magnitude response at extreme frequencies will give you a sense of highpass vs. lowpass for the filter. For example,

$$ |H(e^{i0},e^{i0})| = 2\sqrt{(\cos{0} + \cos{0} - 2)^2} = 2\sqrt{(1+1-2)^2} = 0,$$
and
$$ |H(e^{i\pi},e^{i\pi})| = 2\sqrt{(\cos{\pi} + \cos{\pi} - 2)^2} = 2\sqrt{(-1-1-2)^2} = 8.$$

Of course, the transfer function can be evaluated directly to plot the magnitude response of the filter as well. Here's an example using numpy:

import numpy as np
import pylab as py
from mpl_toolkits.mplot3d import axes3d

def H(z1,z2):
    return z1 + z2 + 1./z1 + 1./z2 - 4.0

n = 100
w1 = w2 = np.linspace(-np.pi,np.pi, n)
mag = np.zeros((n,n))

for i1 in xrange(0,n):
    for i2 in xrange(0,n):
        z1 = np.exp(1j*w1[i1])
        z2 = np.exp(1j*w2[i2])
        mag[i1,i2] = np.abs(H(z1,z2))

fig = py.figure()
ax = fig.add_subplot(111, projection='3d')
X, Y = np.meshgrid(w1,w2)
ax.plot_surface(X, Y, mag, cmap='bone', alpha=.5)
py.show()

Keep in mind that if you use the kernel 2D FFT technique, the resulting magnitudes won't necessarily be centered at zero like the above plot.


what are resources you recommend to learn GPS signal processing?

A good place to start is Van Diggelen and Enge's old course, now on YouTube.
There's a lot of GNSS related information on Wikipedia but I would look at Navipedia first.
Finally, you could try one of open source software receivers, the most comprehensive being gnss-sdr.


Your characteristic polynomial is indeed

$$\lambda^2-4\lambda+4=(\lambda-2)^2$$

If you had two distinct roots $r_1$ and $r_2$ the general solution would have the form you suggested:

$$y(n)=C_1r_1^n+C_2r_2^n$$

However, in your case you have a double root $r=2$ which results in a general solution of the following form:

$$y(n)=C_12^n+C_2n2^n\tag{1}$$

The constants $C_1$ and $C_2$ are chosen to match the initial conditions and the input signal $x(n)$, which in order to compute the impulse response is a delta impulse $\delta(n)$. Assuming $y(-1)=y(-2)=0$ we get

$$y(0)=x(0)-x(-1)=1\\
y(1)=4y(0)+x(1)-x(0)=4+0-1=3$$

Comparing to (1) gives

$$y(0)=C_1=1\\
y(1)=2C_1+2C_2=3\Rightarrow C_2=1/2$$


Signal spectrum when carrier width changes

If the signal is band-limited (at least above some noise level) below half the highest frequency in your initial frequency evaluation set, then you can resample the spectrum to any other frequency spacing by using periodic Sinc (Dirichlet) kernel interpolation.

If there is an integer relationship between the two carrier frequency bin spacings, then you can pre-calculate a polyphase interpolation table using the Sinc function (or a windowed approximation, if that meets your S/N requirements).  Otherwise you can convolve the entire spectrum vector against a computed periodic Sinc or Dirichlet kernel with the appropriate phase offset per computed "carrier" frequency.


Moving average filter for sound in MATLAB

1.) The moving average filter works as a low pass filter. Consequently, it will attenuate the high-frequency part of the signal and passes the low-frequency part of the signal. Therefore, the output signal should sound flat.

2.) If we change b0= 0.5 and b1 = -0.5, the filter works basically as a differentiator, subtracting one sample from another sample. In this case, the filter works as a high pass filter, attenuating the low-frequency part of the signal. Therefore, when listening to the output signal, the high-frequency part of the input signal should appear amplified.


You should realize that when you say "telecommunications" you are talking about a huge field involving lots of different technologies.  Based on the mention of AM/FM modulations I infer that you are referring to digital (though AM/FM are actually analog) wireless transmission techniques (though I may have inferred incorrectly since wired modems can use some pretty sophisticated modulations).  Even restricting ourselves to wireless telecommunications, there are lots of different technologies that we could talk about- modulation types are just one.

Anyway, back to what I think you were asking.

The progression of modulation types went something like the following:

AM (very simple analog)
FM (more sophisticated analog- but uses more bandwidth than AM)
OOK (very simple digital)
FSK (digital version of FM)
PSK
QAM
TCM
DSSS
OFDM

If you google any of the above terms and add "modulation", you should get some useful references.

OFDM is, as far as I am aware, the latest-and-greatest modulation type.  It has the useful feature of being extremely efficient at packing data into a given amount of bandwidth.  Since bandwidth is at a premium for the telecoms, that is very important.  Many of the other modulation types are still used though, and will continue to be used because they have other useful features.


Horizontal and vertical gradients are computed by taking neighbor pixel differences:

$$g_{x}=L(x+1,y)-L(x-1,y)\\g_{y}=L(x,y+1)-L(x,y-1)$$

Gradient magnitude is computed the same way as in your formula:

$$m(x,y)=\sqrt{g_x^2+g_y^{2}}$$

Replacing $g_{x}$ and $g_{y}$ with above will give you the original formula.

Gradients are usually computed by forward or central difference formula, but I think in this case the above is still accurate enough since the image $L$ is smoothed.


Problem to implement a sinusoid oscillator in matlab

OK, since the original code was changed a bit, here is the complete solution (with some useful modifications):

n=1000;
w0=pi/16;
beta=1;
alpha=2*cos(w0);
y=zeros(1,n);
y(1)=0;
y(2)=sin(w0);

for i=3:n
    y(i)=alpha*y(i-1)-beta*y(i-2);
end

stem(y)
% End of code