Calculating the perimeter in binary image

Provided that these are representative images, you can apply a region growing operation, somewhere at the bottom of the image (or with the seed point being user defined) which will help in extracting just the contour you want to measure. It might look something like this:

Then, apply the high pass convolution matrix of:

$$\begin{bmatrix} 1 & 1 & 1 \\ 1& -8 & 1 \\ 1 & 1 & 1 \end{bmatrix}$$

Which will get rid of areas of constant "amplitude" and leave only the sharp transition of the perimeter. The image now might look like:

Now you can either just count the white pixels (there are 2658 in this one) or trace the white curve (again, a form of region growing) by starting from one of the points for which the above convolution kernel provided a high output (that is, one of the points that is exactly on the perimeter) recursively by looking for 8-connected components.

Obviously, tracing the curve will provide a much more accurate estimate because there, you can take into account the $\sqrt{2}$ distances of pixels resting on a diagonal. If you then know your pixel size, you can also convert the estimate from pixel space to actual units of length.

Hope this helps.

(All processing done in GIMP, some white pixels outside of the perimeter where detected as well but can be easily taken out manually (even in ImageJ) or, if tracing the curve would not affect the measurement.)


A good textbook for designing signal filters

There are a lot of books out there, but if you are interested in Control and Signal Processing, I strongly suggest you take a look a Stephen Boyd Lectures from standford:

http://www.youtube.com/watch?v=bf1264iFr-w

There's the first one, the entire course is really valuable and he is a great Teacher.

Appart from That here's a good list of my preferred books on Signal Processing, some are more introductory, and some are more advanced:

Basics:

Signals And Systems - Oppenheim & Willsky
Digital Signal Processing - Proakis & Manolakis
Discrete Time Signal Processing - Oppenheim & Schafer

A bit more Advanced and Related to Control:

Optimum Signal Processing - Sophocles J. Orfanidis
Optimal Filtering - Brian D. O. Anderson & John B. Moore
Optimal Control: Linear Quadratic Methods - Brian D. O. Anderson & John B. Moore

Some on Statistical Processing (Advanced):

Statistical Signal Procfessing - Louis L. Scharf
Linear Estimation - Thomas Kailath, Ali H. Sayed & Babak Hassibi (Watch out this one is really tough)
Spectral Analysis of Signals - Petre Stoica & Randolph Moses

And some about Adaptive processing (Useful in every area of signal processing and control):

Adaptive Filter Theory - Simon Haykin
Adaptive Filters - Ali H. Sayed
Satistical And adaptive Signal Processing - Manolakis, Ingle & Kogon

And of course, the father of all Control Engineering Textbooks:

Linear Systems - Thomas Kailath

From your question i cannot really guess your background, but if you have never studied any topic like LTI systems, Digital Systems or Control, i would say you have a couple of years ahead of you to really learn and dominate the topic. The books i listed are the most referenced books in university courses, and in general the ones that everyone studies. They are for the most part, very theoretical, a lot of math and stuff like that. But it's a must have for reference. If you are looking for more instant satisfaction i would say try to google some cookbooks about digital filters or control, with the basics of PID controllers, linear regulators, etc.

Hope it helps


The short-time Fourier transform is generally a redundant transformation, usually implemented with the same subsampling over every frequency. If the window is well chosen, it is complete: you can invert it and recover any initial signal.

Since it is redundant and complete, it has many perfect inverses. It can be implemented and understood using more generic tools: (oversampled) complex filter banks. Given a window type and length plus the overlap provides you with an analysis filter bank for which you can compute whether it is invertible or not. If it does, you can compute a natural inverse, and optimized inverses too. The overlap-add is just one on the many potential inverses, probably the most common, which often restricts the window choice.

Standard discrete wavelet transforms are filter banks too, with the difference that the subsampling is not the same in each frequency band (or more properly scale). This turns into uneven lengths for each scales. Yet, there exists redundant wavelet implementations that yield "a rectangle array" of coefficients you can work with. The best known schemes are called under different names: shift-invariant or time-invariant wavelets, undecimated wavelets, stationary wavelet transform (SWT), and sometimes cycle-spinning. Its standard reconstruction involves steps similar to overlap-add, except they are more "embedded" due to the different sampling factors over scales. You can use them with any discrete wavelet from a library, or even by designing your own wavelet. The reason is that standard discrete wavelets were designed with non-redundancy in mind, which restricts the choice of wavelets. With redundancy, you wavelet choice increases, as the constraints to fulfill are less stringent. The "ultimate" avatar is the continuous wavelet transform, which admits "almost" every inverse synthesis wavelet. My last sentence is quite lousy, I hope you get the meaning: when a square matrix is invertible, it has only one inverse. When a "rectangular" matrix is left-invertible in a generalized way, it has several inverses.

There seems to be a python implementation of the Stationary Wavelet Transform. You can find a few references in the 2.3.4. Translation invariant wavelets chapter of the linked paper.

It is generally much more robust for detection, denoising or restoration in practical applications (geophysics, non-destructive testing, ultrasounds, vibrations).


use of complex conjugate of weight vectors in beamforming literature

This comes from the fact that the convolution operator for discrete-time linear systems is typically defined as the inner product between the input signal $x_i(k)$ in your example) and the system's impulse response ($w_i$). The definition of the inner product for complex vectors includes the conjugation:

$$
a \cdot b = \sum a_i b_i^*
$$

For the above definition, note that there doesn't seem to be consistent agreement over which of the arguments to the dot operator gets conjugated, so be careful when reading the literature. For the application of discrete-time FIR filters (which your example is), the typical convention is to conjugate the vector of filter taps, hence the $w_i^*$ term in the sum.

Note that this convention is not universal in the FIR filtering literature; I have seen papers where the convolution sum for a system with complex-valued impulse response does not include the conjugate operator. However, I prefer to use the form that you showed in your question.


Equation visualization of complex sinusoid function

The exponent function is directly related to the sine and cosine via Euler equation:

$$e^{i \theta} = \cos(\theta) + i \sin(\theta)$$

For example for function $z=e^{i \pi t}$:

$$z=\cos(\pi t) + i \sin(\pi t) $$

You could plot it in 3D, where x axis is the time, and y & z axes are real and imaginary part of your complex signal, respectively. In 3D that is going to be a helix:

But if you look at projections onto the y=0 and z=0 planes, then you will notice sine:

and cosine waves:

On the other hand, if you look at the x=0 projection, then you will observe a (unit) circle:

You can find the applet with this visualization in Wolfram Demonstration Project.


Deriving the Matrix Inversion Lemma for RLS Equations vs the Woodbury Derivation

It is not clear what are you asking but I will try answer both things.

Deriving the Matrix Inversion Lemma

The Matrix Inversion Lemma goes as:

$$ {\left( A + U C V \right)}^{-1} = {A}^{-1} - {A}^{-1} U {\left( {C}^{-1} + V {A}^{-1} U \right)}^{-1} V {A}^{-1} $$

Deriving it is by utilizing these useful identities:
$$\begin{align}
U + U C V {A}^{-1} U & = U C \left( {C}^{-1} + V {A}^{-1} U \right) = \left(A + U C V \right) A^{-1} U \label{auxEqn001}\tag{1} \\
{\left( A + U C V \right)}^{-1} U C & = {A}^{-1} U {\left( {C}^{-1} + V {A}^{-1} U \right)} ^{-1} \label{auxEqn002}\tag{2}
\end{align}$$

Where $ \ref{auxEqn002} $ comes from $ \ref{auxEqn001} $ by multiplying the 2 last terms by $ {\left(A + U C V \right)}^{-1} $ on the left and $ {\left( {C}^{-1} + V {A}^{-1} U \right)}^{-1} $ on the right.

Then it is straight forward:
$$\begin{align}
{A}^{-1} &= {\left( A + U C V \right)}^{-1} \left( A + U C V \right) {A}^{-1} && \text{As $ {\left( A + U C V \right)}^{-1} \left( A + UCV \right) = I $} \\
& = \left( A + U C V \right)^{-1} \left( I + UCVA^{-1} \right) \\
& = \left( A + U C V \right)^{-1} + \left( A + U C V \right)^{-1} U C V {A}^{-1} && \text{Multiplying the terms} \\
& = \left( A + U C V \right)^{-1} + {A}^{-1} U {\left( {C}^{-1} + V {A}^{-1} U \right)}^{-1} V {A}^{-1} && \text{Utilizing \ref{auxEqn002}} \\
& \Rightarrow {\left( A + U C V \right)}^{-1} = {A}^{-1} - {A}^{-1} U {\left( {C}^{-1} + V {A}^{-1} U \right)}^{-1} VA^{-1}
\end{align}$$

The nice thing is we don't need the Matrix Inversion Lemma (Woodbury Matrix Identity) for the Sequential Form of the Linear Least Squares but we can do with a special case of it called Sherman Morrison Formula:

$$ {\left( A + u {v}^{T} \right)}^{-1} = {A}^{-1} - \frac{ {A}^{-1} u {v}^{T} {A}^{-1} }{ 1 + {v}^{T} {A}^{-1} u } $$

The Sherman Morrison Formula is the MIL where $ C = I $, $ U = u $ and $ V = {v}^{T} $.

Deriving the Sequential form of the Linear Least Squares Estimator

In Sequential Form of the Least Squares Estimator for Linear Least Squares Model I derived the sequential form. To add something new I will extend the derivation to the Weighted Linear Least Squares Estimator.


Different way to separate a particular frequency from a signal

If the frequencies to be rejected are the discrete frequencies given, and the sampling rate can be commensurate with these frequencies, then a moving average filter with the coefficients complex rotated to shift the response to 17 Hz can be done to reject all other integer frequencies as in the plot below.

This is the concept of "heterodyning the coefficients" rather than heterodyning the signal. When we heterodyne a signal, the signal is multiplied by a complex rotator ($e^{-j\omega t})$ which will in the process shift it in frequency according to the rate of rotation. In similar fashion, rotating your coefficients will shift the response of the filter. So you can either move your signal to the filter as in a downconversion followed by a low pass filter... or move your filter to the signal as in this case!
Here is the Matlab/Octave code that implemented this filter.
In this case the sampling rate is 26 Hz with a complex (I and Q) signal path.
z=exp(j*2*pi/26);
coeff=z.^([0:25]*17);
[h,w]=freqz(coeff,26,2048,'whole');
plot(w/pi*13,20*log10(abs(h)))
grid on
axis([0 26 -50 0])

(Note, this is the same frequency response as one bin in the DFT)
The coefficients themselves follow the pattern as shown in the following complex diagram, where the first 5 coefficients are labeled, and the rest can be followed along the trajectory lines shown. Observe that the magnitude of each coefficient is one (as in a moving average filter) but the phase is continuously rotated which causes the frequency response to shift at the rate of rotation.

Update:
As further suggested by @MarcusMüller in the comments, actually doing the signal homodyne followed by a low pass filter may likely be simpler in this case given that the OP has variable conditions for the signal of interest. Assuming the interference to be rejected is still discreet tones with an integer spacing relationship to the signal of interest (the signals need not be integer frequencies, but the multiplies of their spacing must be, as the nulls of a moving average filter will be located at n/T Hz where n is any integer and T is the time duration of the moving average in seconds).
So in this case we would first multiply the signal by $e^{-j2\pi 17 t}$ and then follow the (complex) output with a moving average filter (which means two filters, one for the I (in phase or real) path and one for the Q (quadrature or imaginary) path. In this case the moving average length as a minimum must be:
2Hz,10Hz,17Hz,19Hz and 25Hz
Move 17 to 0:
-12Hz, -7Hz, 0 Hz, 2Hz, 8Hz
And the least common multiple is 1 so a 1 second moving average is required with a commensurate sampling rate. Further as @MarcusMüller has suggested, this moving average structure could be implemented as a CIC decimating filter.  The output would be the "DC Equivalent" of the signal of interest, with the same amplitude and phase as the signal at 17 Hz. Further if the signal that is desired changes in frequency, the phase rotator can be changed in one place rather than updating the taps as I had done.
For further reading see:
Phase rotator: Numerically Controlled Oscillator (NCO) for phasor implementation?
CIC: CIC Cascaded Integrator-Comb spectrum
Note also to be considered as a filtering approach is the 2nd order resonator (or cascade of two of these structures), based on a rotated exponential averaging filter. (see FFT analysis for Vibration Signal). The image of this filter structure is duplicated below (see link for details on what the variables mean and how to compute them):


Sketching the signal always helps in such cases. Let's consider a cosine of period $T_0$, that is, of fundamental frequency equal to $f_0=1/T_0$. That is the signal on the left.

As you can see on the right, $|A\cos(2\pi f_0 t)|$ does not share the same period with $A\cos(2\pi f_0 t)$. The period of the signal you are looking for is $T_0^{\prime}=T_0/2$ and its fundamental frequency is $f_0^{\prime} = 2f_0 = \frac{2}{T_0}$.

Assuming that you're looking for the exponential Fourier Series, the coefficients should be $$X_k = \frac{1}{T_0^{\prime}}\int_{T_0^{\prime}}x(t)e^{-j2\pi kf_0^{\prime}t}dt$$
and for $k=0$, which is usually calculated seperately, it is $$X_0 = \frac{1}{T_0^{\prime}}\int_{T_0^{\prime}}x(t)dt$$

Notice that you can integrate over any time interval of duration equal to $T_0^{\prime}=T_0/2$. I will pick $[-T_0/4, T_0/4]$ since the cosine is always positive in this interval and thus removing the $|\cdot|$ is trivial.

I will help with $X_0$ and you can do the $X_k$s. It is:

$$\begin{align}
X_0 &= \frac{1}{T_0^{\prime}}\int_{T_0^{\prime}}x(t)dt = \frac{2}{T_0}\int_{-T_0/4}^{T_0/4}A\cos(2\pi f_0t) dt \\
&= \frac{2}{T_0}\frac{A}{2\pi f_0}\sin(2\pi f_0 t)\Big]_{-T_0/4}^{T_0/4} = \frac{2A}{2\pi f_0T_0}\Big(\sin\Big(\pi f_0\frac{T_0}{2}\Big) -\sin\Big(-\pi f_0\frac{T_0}{2}\Big)\Big)
\end{align}$$

We know that $T_0 = 1/f_0 \Longrightarrow f_0T_0=1$ and that $\sin(-x)=-\sin(x)$. Simplifying:

$$\begin{align}
X_0 &= \frac{A}{\pi}\Big(\sin\Big(\frac{\pi}{2}\Big) + \sin\Big(\frac{\pi}{2}\Big)\Big) = \frac{2A}{\pi}
\end{align}$$

Now for the $X_k$s, I've already mentioned that you have to use
$$X_k = \frac{1}{T_0^{\prime}}\int_{T_0^{\prime}}x(t)e^{-j2\pi kf_0^{\prime}t}dt$$
You can use the same interval for your integration, keep in mind that $f_0T_0=1$ to simplify your expressions, and treat $k$ as a constant until the end. Judging from your question, you should get something that includes $\frac{2A}{(2k+1)\pi}$.


What are the options for turning a bit stream into a byte stream?

I would use a frame alignment word to synchronize your data.  The basic idea is to insert your FAW at the beginning of every $N$ bytes.  After you demodulate the data you look for the FAW to synchronize the data.  To avoid false positives you should try to have a longish FAW (16 bits should be plenty) and look for repeated FAWs.  To avoid inefficient overhead you should try to make $N$ large compared to the size of the FAW.


Again I refer you to the Murray's book. Which is free by the way

http://www.cds.caltech.edu/~murray/books/AM08/pdf/fbs-limits_18Aug2019.pdf

From the book "Furthermore a zero is said to be “slow” if its magnitude is smaller
than the intended closed loop bandwidth."


I'm wondering after finding the magnitude squared values at each FFT bin, do you then add each bin value up and divide by the number of bins (N) to get the average power?

Yes.

  And also, do you just repeat this averaging M times?

Yes, you would have to...

  And why exactly do you repeat this averaging M times??

...to reduce the impact of noise on your measurement and also, depending on the length of the FFT, to make sure that you have captured the whole dynamics of your signal.

For more information, please see Power Spectrum Estimation Using the FFT.

Hope this helps.


MATLAB, How to filter a discrete signal?

Easiest way to remove baseline is to remove the average:

filtered = original - mean(original);

Actually, the average is the first coefficient of the Fourier transform, so it is actually a very sharp filtering: you are eliminating the DC coefficient.

If you want more classic filtering, then check functions like butter and friends, which will synthesize an IIR filter, then use filter to filter out your signal.

Matlab also includes a filter design tool.


How to index output of convolution product

If $x[n]$ and $y[n]$ are both causal and starting at index $0$, then the result of convolution will also be causal and it will start at index $0$. Just plug in $n=-4$ in the expression for $z[n]$, you will find that it will be $0$.
$$z[n] = \sum^{\infty}_{k=-\infty}x[k]y[n-k]$$
First non-zero term in above expression is at $n=0$. Because for $k<0$, $x[k]$ will be $0$, you can re-write the above expression as:
$$z[n] = \sum^{\infty}_{k=0}x[k]y[n-k]$$
In the above expression, put $n<0$, you will see that $y[n-k]$ will always be 0.
There will be a transient delay of $N-1$ samples, when you are using Moving-Average FIR Filter of length $N$, not an advance. The delay is due to the fact that you need $N-1$ previous samples of the input before you can produce an average of $N$ samples. You cannot have an output before you feed input to a causal FIR filter. Expecting an output at $-21sec$, is expecting to look into future which is not correct here.
The figures you have posted is showing that you have 8399 samples of output starting at $n=0$. As you can see, approximately at $5sec$, the seismic activity starts, Moving average just smoothens-out the original input signal with a transience of $0.3sec$.


As mentioned already, it is incredibly easy to read metadata off of WAV files as it is a tag based format and there are already implementations of readers for this information out there (e.g. tinytag in Python).

Having said this though, "bpm, key, instrument" are not part of typical information you would find in a WAV file. If the creator of the file put them there, you will be able to read it, but "bpm" is not essential for playing back an audio file. What you are more likely to be able to retrieve is things like Sampling Frequency, Bits per sample, number of channels and of course anything to do with any compression method used. Building an automatic "recogniser" which could be populating these fields given some WAV file, is possible but not exactly straightforward.

Bpm, key and instrument however are essential for playing back MIDI files. Perhaps these are more useful to you (?)

Hope this helps.


How to compute impulse and frequency response of Flanger?

Only linear time invariant systems are fully characterised by their impulse response. Modulation effects like a flanger break time invariance. While you can still create a momentary impulse response (which is only of very limited use), there is no meaningful frequency response.


CIC SINC3 Filter with unsigned Integrators

Addition and subtraction in 2's complement or unsigned is the same. The only difference is for overflows.

Assuming 16-bit arithmetic, overflows occur in unsigned addition when you go over 0xFFFF (65535). For example, if you add 1 to 0xFFFF you will wrap around to 0.

In 2's complement, the overflow will technically occur when you go over 0x7FFF (32767). If you add 1 to 0x7FFF you will wrap around to 0x8000 (-32768).

Notice however that 0xFFFF + 0x0001 = 0x0000 and 0x7FFF + 0x001 = 0x8000 whether you're in 2's complement or unsigned.

Therefore, it does not matter whether it is implemented in 2's complement or unsigned, the behaviour will be the same.


Discrete Fourier transform - finding the fundamental quickly?

What you are describing is very similar the the Harmonic Product Spectrum method of pitch estimation, as listed in this Stanford CCRMA paper.

An FFT does not give you an "infinite sum of amplitudes", but a finite number of result bins depending on the length of the FFT.

5 mS is only 1 period of a 200 Hz note, and only a fraction of a period below 200 Hz.  Musical pitch recognition usually requires hearing or analyzing a multiple number of periods of the periodicity of a pitched sound.  And a lot of music uses notes below G2. If you have a sufficient length of data, computing a pitch estimate from that data might take only on the order of microseconds rather than milliseconds on a modern PC or mobile device.


It seems I found an answer on the Panotools Wiki:

  Unlike rectilinear lenses, fish-eye lenses do not follow the
  tangent-plane geometry, but instead have built-in distortions designed
  to achieve wide fields of view. The radial lens distortion parameters
  are used the same way for rectilinear lenses and fisheye lenses, but
  they should never be used to attempt to remap a fisheye to a
  rectilinear image. This is done by selecting the proper source and
  destination projection. Fisheye geometry follows a rapidly-changing
  trigonometric function which can hardly be approximated by a third
  degree polynomial.

  For fisheyes, the lens correction parameters correct for the deviation
  between a real lens and the ideal fisheye geometry.

The last sentence suggests that the lens correction is done to get a proper geometry of whatever lens it is. The fisheye is therefore remapped after the lens corrections are done.

Looking at the Fulla tool source code from Hugin, it seems that it performs corrections in the following order:

Flat-field correction (i.e. removing vignetting
Radial distortion correction (barrel, pincushion, wavy distortion)
TCA correction (i.e. removing one type of chromatic aberration)


On of the easiest ways to see why the outcome of the problem that Wiener Filter tries to solve is Linear Time Invariant (LTI) Filter is to derive (Solve) it.

By the very definition of the Wiener filter we're after a solution in the form of a filter applied using convolution which means this filter must be LTI.

Let's go through this.
The Wiener Deconvolution problem is given by:

$$ y \left( t \right) = \left( h \ast x \right) \left( t \right) + n \left( t \right) $$

Where:

$ x \left( t \right) $ is the signal to be estimated.
$ h \left( t \right) $ is the known impulse response of a linear time invariant system.
$ n \left( t \right) $ is some unknown additive noise, independent of $ x \left( t \right) $.
$ y \left( t \right) $ is our observed signal.

The Wiener deconvolution, in order to solve the above problem (Estimating $ x \left( t \right) $) is finding $ g \left( t \right) $ such that:

$$ \hat{x} \left( t \right) = \left( g \ast y \right) \left( t \right) $$

Minimizes:

$$ \mathbb{E} \left[ { \left( x \left( t \right) - \hat{x} \left( t \right) \right)}^{2} \right] $$

Since we defined the Winer Filter $ g \left( t \right) $ by the Convolution Operation we defined it as LTI.


Have a look at this plot from this page at the Mathworks site:

The main issue is, as that page says:

  The frequency input to the Chebyshev Type II design function sets the beginning of the stopband rather than the end of the passband.

As you can see, the Chebyshev Type I filter has controllable passband gain out to the cutoff frequency.  That is not true for the Type II: Only the stopband gain is controllable there.  The passband roll-off for Type II makes it not very interesting for most purposes.


Phase shift of all-pass bi-quadratic filter

First of all, that formula for the phase has a sign error, at least if $\alpha$ is defined as in the diagram of the filter structure. Second, your phase plot shows a jump from $-\pi$ to $\pi$, which is no jump at all because a phase of $\pi$ is the same as a phase of $-\pi$. The jump occurs only because the arctangent function returns the principal value of the phase in the interval $[-\pi,\pi]$. You can get rid of the jump by using the phase unwrapping function unwrap in Matlab/Octave. But note that this is just a cosmetic thing, there is no phase jump in the physical sense of the word.

Now for the derivation of the formula for the phase. From the given filter structure we get the following equation for the $\mathcal{Z}$-transforms of the input and output signals:

$$Y(z)=X(z)z^{-2}+\alpha\left(X(z)-Y(z)z^{-2}\right)\tag{1}$$

From this the transfer function can be computed as

$$H(z)=\frac{Y(z)}{X(z)}=\frac{\alpha+z^{-2}}{1+\alpha z^{-2}}=\frac{\alpha z+z^{-1}}{z+\alpha z^{-1}}=\frac{B(z)}{B\left(\frac{1}{z}\right)}\tag{2}$$

with $B(z)=\alpha z+z^{-1}$.
$H(z)$ as given by $(2)$ is indeed an allpass transfer function with $|H(z)|=1$ for $|z|=1$. Note that the system described by $(2)$ is only stable for $-1<\alpha <1$.

On the unit circle $z=e^{j\omega}$, $H(z)$ can be written as

$$H(e^{j\omega})=\frac{B(e^{j\omega})}{B^*(e^{j\omega})}\tag{3}$$

where $^*$ denotes complex conjugation. Consequently, the phase of $H(e^{j\omega})$ is given by

$$\phi_H(\omega)=2\phi_B(\omega)\tag{4}$$

where $\phi_B(\omega)$ is the phase of the numerator $B(e^{j\omega})$:

$$\begin{align}\phi_B(\omega)&=\arg\big\{\alpha e^{j\omega}+e^{-j\omega}\big\}\\&=\arg\big\{(1+\alpha)\cos(\omega)-j(1-\alpha)\sin(\omega)\big\}\\&=-\arctan\left(\frac{1-\alpha}{1+\alpha}\tan(\omega)\right)\end{align}\tag{5}$$

The phase of $H(e^{j\omega})$ is thus given by

$$\phi_H(\omega)=-2\arctan\left(\frac{1-\alpha}{1+\alpha}\tan(\omega)\right)\tag{6}$$

which differs from the given formula in the sign of $\alpha$.


Is R suitable for digital signal processing

Since the bulk of R’s DSP capability comes from the signal package which was ported over from the open source project Octave (itself influenced by MATLAB), there's no intrinsic limitation of R.

What you have picked up on, are ecosystem preferences. We learned MATLAB in university, picked up numpy/scipy/sklearn at work, so R isn't the first weapon of choice.

That said, there are some interesting "bridge" projects that you may be interested in exploring: https://stackoverflow.com/questions/5630441/how-do-rpy2-pyrserve-and-pyper-compare.

There is also the excellent interactive python (iPython) project, now called Jupyter, that recently introducing support for R notebooks: https://www.continuum.io/blog/developer/jupyter-and-conda-r. So tools & their inter-operability are converging!


Who first understood the importance of poles?

If you consider poles of an integral transform domain to be important to the solution of differential equations: (as usual,) Euler did it first, 1753.
One "importance" of poles is that they're part of a very useful representation for linear systems. They must've appeared when people started looking at functions as built from generating functions, so that'd date that .
So, both Bode and Nyquist are probably 100 years too late: this would emerge from looking at systems like Laplace did, around 1785.
If you think of poles as something that is relevant to your analysis of holomorphic functions, then clearly Cauchy would be your guy, starting with "Sur les intégrales définies", 1814, where he looked at simple poles within an area of integration.
When you really look into where the residue theorem (which really concerns itself with integrals over functions with poles, and is immensely important to a lot of things that we can do with the Fourier transform) was formed, you'd be looking at Lindelöf's "Le calcul des résidus et ses applications à la théorie des fonctions", 1905.
Finally, the Nyquist stability criterion was published when he was working at Siemens, 1930.

     Laplace     Lindelöf
  Euler | Cauchy    |
    |   |  |        |  Nyquist
    v   v  v        v  v
|---------.---------.---------.-->
1         1         1         2
7         8         9         0
0         0         0         0
0         0         0         0


Acquisition of vibration signal and processing help

I suppose that you are interested in measuring vibration for monitoring purposes, whether for the shaft's bearings or what the shaft itself is doing as it gets loaded (e.g. sagging and causing vibration).

For this purpose, the vibration signal can be assumed stationary. This means that given enough time, no new components will be added to the vibration signal. The vibration is "the same" in every rotation.

With these assumptions in place, what you can do is to integrate an aspect of your vibration signal over short periods of time and cross referencing it with your RPM.

For example, let's say that the vibration signal that you get from your sensor is represented by some $x(n)$. If you want an estimate of how strong the vibrations get, you can simply integrate the $|x(n)|$ over time (that is, the absolute value of $x$) and see if the amplitude of the vibration seems to increase with respect to RPM. The integration itself can be done with a simple moving average filter. If your Sampling Frequency is 7kHz ($F_s > \frac{f_h}{2}$) and your highest frequency $f_h$ here seems to be 2.5kHz), then this means that 1 revolution at 250Hz takes about 140 samples or approximately 20ms. So, if you were to filter your $|x(n)|$ with a simple moving average filter with length 140 samples, you would generate another signal that would give you average amplitude of vibrations per revolution at your low end RPM. As the RPM increases of course, the smoothing will be greater.

But, you are also very lucky because in addition to your vibration signal you also have the RPM sensor. From these two and the point of view of monitoring the bearings you can understand two things: Whether the bearing's shell has developed a chip, or whether one of the bearing elements is deformed. Here is why: If the bearing shell has a chip then there will be a peak in the vibration signal at the same phase in every revolution. If one of the bearing elements is defored (whether cylinder, ball, other) then, as the bearing rotates, the elements meet the shell at a different point at each revolution. Therefore, your get a "rolling" pulse in the signal. On the first revolution it occurs 2ms after the rotation starts, on the second revolution it occurs 4ms after the rotation starts, on the third revolution it occurs 6ms after the rotation starts and so on.

To see these things, you can plot your vibration signal in parallel (or, in phase) with your rotation signal and look for "regular bumps" or "rolling bumps".

The same applies for when you want to measure if your axle is sagging beyond an RPM (because it was asymetrically lathed or it was loaded suddenly or whatever). You just monitor the vibration signal for bumps versus the RPM.

Towards that direction, you might find wavelet analysis more useful.

Hope this helps.


How to save Fourier series approximated signal to a WAV file

OK, after digging the case by myself I managed to generate square wave files for both cases. Updated the source code found in my opening post.
In case of combined Fourier Series and Fejér taper approximations of square wave signal, result looks like this (plot done with Audacity):

and spectrum like this (Plot done with Octave):

(green = w_sqr, black dotted = Ffej, fc=1230Hz, fs=44100Hz, w_sqr and Ffej are data holders from source code)


Audio processing -- difference between block and sample-by-sample?

Note that FIR filters are sometimes implemented using an FFT (overlap-add, overlap-save). In that case it makes sense to have power of 2 buffers lengths (depending on the FFT implementation). This is of course an example of block processing, where you have to wait for a whole block before the computation of the output signal can begin. The consequence is that you always have some latency. Sample by sample processing is possible with a time-domain implementation, where you get one output sample for each input sample. Note that a time-domain implementation is not necessarily sample-by-sample but could also be using block processing.


Non Periodic, Deterministic Power Signals

Consider the signal
$$x(t) = \sum_{n=-\infty}^\infty \operatorname{rect}(t-n)\sin(2\pi nt)$$ which consists of $|n|$ periods of the signal $\sin(2\pi nt)$ in the interval $\left(n-\frac 12,n+\frac 12\right)$ for each integer $n$. Clearly, $x(t)$ is not periodic. Equally clearly, $x(t)$ is a power signal with average power $\frac 12$, and thus it is a nonperiodic but deterministic power signal. I leave it to the OP to determine whether the signal in question has a Fourier transform in the generalized sense (meaning it involves impulses) or not.


Consider one-dimensional function $f(x)$.

The first order taylor expansion is $f(x_0+h) \approx f(x_0) + f'(x_0)h$

The second order taylor exapnsion is $f(x_0+h) \approx f(x_0) + f'(x_0)h + \frac 1 2 f''(x_0)h^2$

Now we expand three-dimensional function.
$$ D(\mathbf{x_0}+\mathbf h) \approx D(\mathbf{X_0}) + \bigg(\frac{ \partial D}{\partial \mathbf x}\bigg)^T\bigg|_{\mathbf x=\mathbf x_0}\mathbf h + \frac 1 2 \mathbf h^TH(\mathbf x)\mathbf h $$
, where
$$
\frac {\partial D}{\partial \mathbf x} = \begin{bmatrix} \frac {\partial D}{\partial x} \\ \frac {\partial D}{\partial y} \\ \frac {\partial D}{\partial \sigma}\end{bmatrix} =
\begin{bmatrix} \frac {D(x+1,y,\sigma) - D(x-1,y,\sigma)}{2} \\ \frac {D(x,y+1,\sigma)-D(x,y-1,\sigma}{2} \\ \frac {D(x,y,\sigma+1) - D(x,y,\sigma-1)}{2} \end{bmatrix}
$$
$$ H(\mathbf x) = \begin{bmatrix} D_{xx} & D_{xy} & D_{x\sigma} \\
D_{yx} & D_{yy} & D_{y\sigma} \\
D_{\sigma x} & D_{\sigma y} & D_{\sigma\sigma} \end{bmatrix} $$

Additionally, $D_{xy}$
$$= \frac { \frac{D(x+1,y+1,s)-D(x-1,y+1,s)}{2} - \frac{D(x+1,y-1,s)-D(x-1,y-1,s)}{2} } 2 $$

Let's go back.
$$D' = \bigg ( \frac{ \partial D}{ \partial \mathbf x} \bigg ) ^T + H(\mathbf x) \mathbf h = 0 $$ $$\mathbf h = - H^{-1}(\mathbf x) \bigg ( \frac{ \partial D }{\partial \mathbf x} \bigg ) ^T $$

Assume that $\mathbf{x_0}$ is discrete feature point. Iteratively we can get $\mathbf h$.

If $\mathbf{x_0}$ is the closest point to $\mathbf h$, then we select $\mathbf{x_0}+\mathbf h$. (Imagine that 3-dimensional space)

My english ability is poor and there may be some errors.


Video stabilization using Lucas-Kanade method

For video stabilization, you need to estimate motion vector with respect to a reference frame $F_{reference}$ (which let's assume it to be the first frame). Then your motion model should be determined, for simplicity let's assume it to be Translational motion. Therefore, we need to estimate motion in $x$ direction and motion in $y$ direction for each frame. Now you have two options to estimate the change. If you are interested to estimate it through optical flow you can do either these:

Estimate regular conventional optical flow from frame to frame
Estimate Optical flow of first frame to the $i$ $th$ frame.

  I accumulate during the whole process- and practically...

If you've chose the first, then Yes, you need to accumulate the optical flow, but I think it would be simpler to go for the second method.

As a hint, I suggest track a token or compute optical flow of a small patch of the image.

Edit:

I feel the problem is with your Warp function, the code below is part of video stabilization assignment in my class, try something like this sample. The translational parameters are estimated using token based tracking.

function [ ou_I ] = frame_translation_func( in_I,Tx,Ty,ou_RoWs,ou_CuLoMns,ou_layer)
ou_I = uint8(zeros(ou_RoWs,ou_CuLoMns,ou_layer));
for i2=1:ou_RoWs
    for j2=1:ou_CuLoMns
        if (i2 - Tx<1) || (j2 - Ty<1) || (i2 - Tx>ou_RoWs) || (j2 - Ty>ou_CuLoMns)
            ou_I(i2,j2,:)=0;
        else
            ou_I(i2,j2,:) = in_I(i2 - Tx,j2 - Ty,:);
        end
    end
end
end


time-domain channel estimation based on two vectors optimization

You specify that $T$ is known and invertible, so you know $X_1, X_2$, and then it's really trivial: $y[3]=[X_1,X_2]*[h_1, h_2]$, and $y[9]=[X_7, X_8] *[h_1,h_2]$; write that down as matrix system
$$\begin{pmatrix}y[3]\\y[9]\end{pmatrix}=\begin{pmatrix}X_1&X_2\\X_7&X_8\end{pmatrix}\begin{pmatrix}h_2\\h_1\end{pmatrix}$$
and solve it.


Prove that a specific auto-correlation function is even

Let me get into more additional detail from what Peter wrote.
First, note that $h[n]$ is the output of an IIR filter with input $G\delta[n]$. Now, note that such system is causal (uses only past information), and that the input is the impulse function. We know that if the system is causal, the impulse response of such system is causal too. Thus, $h[n]$ is a causal signal: $h[n]=0,\forall n<0$.
Now, let's get into the autocorrelation function thing. From your procedure you have the following:
\begin{align}
\hat{R}(m)&=\sum_{n=m}^{\infty}h[n-m]h[n]\\
\end{align}
So, you want to write this summation starting from $n=0$ right? If $m=0$ we are done. There are other two case: $m<0$ and $0<m$. In your procedure you did:
\begin{align}
\sum_{n=m}^{\infty}h[n-m]h[n] = \sum_{n=0}^{\infty}h[n-m]h[n] - \sum_{n=0}^{m-1}h[n-m]h[n]\\
\end{align}
However, this basically assumed that $0< m$. Otherwise, the second summation does't make sense since the lower limit would be greater than the upper limit. What if $m<0$? Then you would have:
$$
\sum_{n=m}^{\infty}h[n-m]h[n] = \sum_{n=0}^{\infty}h[n-m]h[n] + \sum_{n=m}^{-1}h[n-m]h[n]
$$
Lets inspect both cases under the causality property of $h[n]$. Hence, for the first case (with $0< m$) we have:
$$
\sum_{n=0}^{m-1}h[n-m]h[n] = 0, \ \ \text{since } h[n-m]=0 \text{ for } n-m<0
$$
this is true, since from the summation index range we have $0\leq n \leq m-1$ or equivalently $n<m\implies n-m<0$ (thus $h[n-m]=0$).
Now, for the second case (with $m<0$) we have:
$$
\sum_{n=m}^{-1}h[n-m]h[n]=0 , \ \ \text{since } h[n]=0 \text{ for } n<0
$$
this is true, since from the summation index range we have $m\leq n \leq -1$ or equivalently $n<0$.
Thus, no matter what value of $m$ you have, the result is
$$
\hat{R}(m)=\sum_{n=0}^{\infty}h[n-m]h[n] = \hat{R}(-m)
$$
Note that this works due to the causality property of $h[n]$: this same reasoning works for any causal signal, regardless of its particular definition (as the one you gave with $h[n]$).


Speech Activity Detection- Adhoc Problem with multi source sounds(birds,trees,heavy air) and how to eliminate them completely

This is not aliasing this is not any phenomenon, there are number of birds whose vocal fall in speech range even a frog vocal falls in speech range so essentially a training based speech detection is required to classify whether a given frame is speech or a non speech

one such implementation is shown here

Speech Discrimination based on multiscale spectro temporal modulations


As with all ocr, you need to segment each individual character. Assuming this step is already done (and maybe some post processing to One easy way is to take the number and reduce the resolution. Similar to a lcd character display.

You can then do a simple correlation between the reduced resolution number from the image, and a template. Since the numbers on street signs are printed and very consistent a simple template match could work well.

your template should contain the numbers 0-9 in a font that closely matches those used on your road signs.


How to calibrate a signal with respect to frequency

Potentially a lot of problems.

You really want to apply a filter to your data that the calibration curve represents, which seems to be your conceptual approach, only frequency domain filtering doesn't involve overlapping and windows, which one does in spectral estimation.

The hard part is to convert your calibration curve into a filter.


Clarification on the application of window pre-FFT

Generally, any window (rectangular or Von Hann) convolves the input with the transform of the window.

The transform of a rectangular window is a Sinc shaped "ringy" waveform with infinite ripples, but these side lobes are invisible for pure sinusoids that are exactly integer periodic in the FFT width.

Smoother windows have lower ripples, but in exchange have wider main lobes.


removing vertical lines in grayscale images by thickness

Here is some Python code that does mostly what you want.  You need to add the thick/thin line selection logic.  You'll probably want to tweak it some more too.  There seem to be some artifacts when it is turned back into a .jpg.

from PIL    import Image

import numpy as np

#============================================================
def main():

#---- Load the Image

        theImage = Image.open( "test.jpg")

        thePixels = theImage.load()

        theRangeMin = 0
        theRangeMax = 7

#---- Get the Size

        theSize = theImage.size

        theWidth  = theSize[0]
        theHeight = theSize[1]

        print theWidth, theHeight

#---- Scan some Rows

        theSum = np.zeros( theWidth )
        theScanCount = 0

        for n in range( 100, theHeight, 100 ):
            print n
            theScanCount += 1

            for m in range( theWidth ):
                theShade = thePixels[m,n]
                theSum[m] += theShade[0] + theShade[1] + theShade[2]

#---- Calculate Derivatives

        theFirst  = np.zeros( theWidth )
        theSecond = np.zeros( theWidth )
        theSum /= ( 3.0 * theScanCount )

        for m in range( 1, theWidth-1 ):
            theFirst[m]  = ( theSum[m+1] - theSum[m-1] ) * 0.5
            theSecond[m] = theSum[m+1] - 2.0 * theSum[m] + theSum[m-1]
            print m, theSum[m], theSecond[m], theFirst[m]

#---- Find Vertical Lines

        theThreshold = 10

        theLastFirst = 0

        theLines = []

        m = 1
        while m < theWidth-1:
            if theSum[m] < theThreshold:
                while m < theWidth-1:
                    if theSum[m] > theThreshold and theFirst[m] < theThreshold:
                        break
                    m += 1

                print theLastFirst, m
                theRange = m - theLastFirst
                theRate = 1.0 / float( theRange )
                if theRangeMin <= theRange and theRange <= theRangeMax :
                    theLines.append( [ theLastFirst, m, theRange, theRate ] )

            elif abs( theFirst[m] ) < theThreshold : theLastFirst = m

            m += 1

#---- Hide the Lines

        theLineCount = len( theLines )

        for n in range( theHeight ):
            for p in range( theLineCount ):
                theLine = theLines[p]

                theSpot = theLine[0]
                theFirstShade = thePixels[theSpot,n]
                theLastShade  = thePixels[theLine[1],n]
                theRange = theLine[2]
                theRate  = theLine[3]
                theFirstVal = ( theFirstShade[0] + theFirstShade[1] + theFirstShade[2] ) / 3.0
                theLastVal  = ( theLastShade[0]  + theLastShade[1]  + theLastShade[2] ) / 3.0
                theSlope = ( theLastVal - theFirstVal ) * theRate

                if n == 100: print theSpot, theFirstVal, theLastVal

                for r in range( 1, theRange ):
                    theVal = theFirstVal + r * theSlope
                    theGrey = int( theVal ) & 0xFF
                    theShade = ( theGrey << 16 ) + ( theGrey << 8 ) + theGrey
                    thePixels[theSpot+r,n] = theShade
                    if n == 100: print "     ", theVal

#---- Write the Image to a File

        theImage.save( "out.jpg" )

#============================================================
main()

Followup:

Here is the debug output at the first thin line:

194 119.857142857 8.85714285714 -1.14285714286
195 123.142857143 -3.71428571429 1.42857142857
196 122.714285714 -64.7142857143 -32.7857142857
197 57.5714285714 13.8095238095 -58.2380952381
198 6.2380952381 88.0 -7.33333333333
199 42.9047619048 42.7142857143 58.0238095238
200 122.285714286 -79.2380952381 39.7619047619
201 122.428571429 0.428571428571 0.357142857143
202 123.0 -0.142857142857 0.5
203 123.428571429 0.0 0.428571428571
204 123.857142857 0.142857142857 0.5
205 124.428571429 -0.285714285714 0.428571428571
206 124.714285714 1.42108547152e-14 0.285714285714


Calculating energy spectral density function of half-cosine function

For a "deterministic finite energy" continuous time signal $x\left(t\right)$, as you have exemplified, the energy density spectrum, or equivalently energy spectral density, which is a real nonnegative signal, meant to be used as a function showing the energy distribution of the signal with respect to the frequency is defined to be: $$S_x\left(\omega\right) = \left|X\left(j\omega\right)\right|^2 $$  where $X\left(j\omega\right)$ is the Continuous-Time Fourier Transform $\left(CTFT\right)$ of the signal $x\left(t\right)$ given as: $$X\left(j\omega\right) = \int_{-\infty}^{\infty}{x\left(t\right)e^{-j\omega t} dt}$$

Now coming to the specific sample of your question, you should first find the CTFT $G\left(j\omega\right)$ of $g\left(t\right)$: That your signal is a multiplication of two other signals you should use a fundamental property of CTFT namely the ${multiplication}$ in time equals ${convolution}$ in frequency. Hence

$$x\left(t\right)y\left(t\right) \leftrightarrow \frac{1}{2\pi} X\left(j\omega\right) \star Y\left(j\omega\right) $$

From a table of FTs (or by carrying out the integral if you can) find that for
$x\left(t\right)= \textrm {rect} \left( \frac {t}{T} \right)$ which is defined from $t=-T/2$ to $t=T/2$, we have: $$X\left(j\omega\right) = \frac{2 \sin \left({\omega \frac{T}{2} }\right)} {\omega}$$ and for the $y\left(t\right) = \cos\left( \pi \frac {t}{T}\right)$ we have $$Y\left(j\omega\right) = \pi \delta \left(\omega - \frac {\pi}{T}\right) + \pi \delta \left(\omega + \frac{\pi}{T}\right)$$

What remains now is to compute their convolution, which is easy considering the convolution with impulse given as:$X\left(j\omega\right) \star \delta\left(\omega - \omega_0\right) = X\left(j\left(\omega - \omega_0\right)\right) $

Hence we get: $$G\left(j\omega\right) = \frac{1}{2\pi} X\left(j\omega\right) \star Y\left(j\omega\right)$$

as $$\frac{\sin \left({\left(\omega- \frac {\pi}{T}\right) \frac{T}{2} }\right)} {\left(\omega- \frac {\pi}{T}\right)} + \frac{\sin \left({\left(\omega + \frac {\pi}{T}\right) \frac{T}{2} }\right)} {\left(\omega + \frac {\pi}{T}\right)}$$

which, after simplyfying, becomes:
$$G\left(j\omega\right) =  \cos\left({\omega\frac{T}{2} }\right)  \left( \frac{-1} {\left(\omega- \frac {\pi}{T}\right)} + \frac{1} {\left(\omega + \frac {\pi}{T}\right)} \right)$$

whose magnitude is $$\left|G\left(j\omega\right)\right| = \left|\cos\left(\omega \frac{T}{2}\right)\right| \cdot \left|\left( \frac{ \frac{-2\pi}{T}}{\omega^2 - \frac{\pi^2}{T^2}} \right)\right|$$

And from which we conclude for the Energy Spectrum Density to be:
$$S_x\left(\omega\right) = \left|G\left(j\omega\right)\right|^2 = A^2 \cos^2\left(\omega \frac{T}{2}\right)  \left( \frac{ \frac{2\pi}{T}}{\omega^2 - \frac{\pi^2}{T^2}} \right)^2 $$

where I have added the last linear scaler A in $g\left(t\right)$.


The very first thing to try would be to implement a rough speech/music classifier to detect a radio host/DJ; and to extract the amplitude envelope to detect fade outs and fade ins.

Once you have implemented that and/or if this doesn't solve your problem because your stream is continuous (DJ-style mix), you'll certainly run into the song boundary / song section boundary ambiguity - but if this is not a problem to you then there are solutions. What you are looking for is known in the academic community as ``unsupervised audio stream segmentation" - with a variety of techniques proposed over the years. You have to decide on two ingredients: which audio features to extract from the audio stream; and which statistical analysis method to use to detect changes/structure.

Regarding the features, the consensus is on a combinations of the following:

MFCC (which carry timbre/instrumentation information)
Chroma (which carry key/chord information)
Loudness
Tempo or other measures of rhythm (beat histogram, Pampalk's fluctuation patterns, etc.)

Very frequently, an algorithm can be repurposed just by a change of features... For example, there are papers on musicological analysis of piano music which will not use MFCC (same timbre throughout the piece), or papers on pop song structure analysis which will not use tempo (assumed to be constant throughout the piece, so not informative) - but the methods they describe could very well work for you if you decide to include other features - so don't be narrow when looking for literature.

Regarding the statistical analysis methods, the usual suspects are:

Computation of a self-similarity matrix, and "checkerboard" filter to detect blocks. The classic paper is by Foote, "Automatic Audio Segmentation Using A Measure of Audio Novelty". Foote's novelty detection function can be computed online and incrementally - over a buffer of past audio frames.
Hidden Markov models, used in an unsupervised way. You hope that the HMM states will map well into segments (Peeters and Rodet, "Music Structure Discovering Using Dynamic Audio Features for Audio Summary Generation: Sequence and State Approach"), or you overestimate the number of HMM states and look for changes in the HMM state distribution within small windows (Levy and Sandler, "Structural Segmentation of Musical Audio by Constrained Clustering"). These methods are more geared towards offline processing, and situations in which you have an idea in advance of how many segments you want to identify.
"Streaming" change detection methods. The principle is to split the past $N$ frames of audio into two $N/2$ frames blocks $S_1$ and $S_2$, and to perform a statistical test to decide which of the hypotheses "Frames from $S_1$ and $S_2$ are drawn from the same distribution $P$" and "Frames from $S_1$ and $S_2$ are drawn from a different distribution $P_1$ and $P_2$" is the most likely. Big peaks or dips in the test outcome indicate that the middle of the block is a segment change. There are a number of information criteria used to weigh the two hypotheses, such as the BIC (Classic paper: "Unsupervised Audio Stream Segmentation And Clustering Via The Bayesian Information Criterion" by Zhou). The limitation here is that you usually have little data to perform the test, and the data has weird-looking distributions. See the section 2.2 of Gillet et al's "Comparing audio and video segmentations for music videos indexing" for a few pointers on kernel methods that overcome these limitations.

For a real-life application (indexing of radio streams), see for example Ramona's "Combined supervised and unsupervised approaches for automatic segmentation of radiophonic audio streams" - this work combined the kind of unsupervised section detectors discussed here to a supervised speech/music classifier (to tackle the most obvious part of the problem).


estimate the impulse response after equalization

Note that what you're trying to do is equalization, as opposed to channel estimation. If we ignore the noise for the moment then, ideally, the concatenation of the equalizer and the channel (modeled as a linear system) should be a pure delay:

$$(h\star w)[n]\stackrel{!}{=}\delta[n-K]\tag{1}$$

where $h[n]$ is the channel impulse response, $w[n]$ are the coefficients of the linear equalizer, and $K$ is some delay that depends on the system design. So by using the LMS algorithm to adapt a linear equalizer you do not directly estimate the channel coefficients $h[n]$. As far as I know there is no straightforward way to reliably estimate the channel weights directly from the equalizer weights.

Regarding your second question, the LMS algorithm doesn't require the source signal to be Gaussian.

A final comment on your code: with only $2$ filter taps you won't be able to sufficiently equalize the given channel. And I would recommend to use vectors in the coefficient update, so you can use different filter lengths without changing your code.


Discrete Fourier Transform frequency bin size

The sum in the definition of a length $N$ DFT always goes from $n=0$ to $n=N-1$:

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

However, if the sequence $x[n]$ is real-valued, which is the case for most applications, then the $N$ DFT bins $X[k]$ are not independent of each other:

$$X[k]=X^*[N-k],\qquad k=0,1,\dots,N-1\tag{2}$$

So for even $N$, only the first $N/2+1$ bins carry information, the remaining $N/2-1$ bins can be computed from $(2)$.

Note that for real-valued $x[n]$, the values $X[0]$ and $X[N/2]$ are real-valued, and the remaining $N/2-1$ values for $k=1,2,\ldots,N/2-1$ are generally complex-valued. This means that $N$ real-valued samples $x[n]$ are represented by $N$ real-valued numbers in the frequency domain (2 real numbers plus $N/2-1$ complex numbers). This shouldn't come as a surprise.


This is indeed correct. There is a few things you could do to make it a better graph:

Label the X-axis
Label the Y-axis
Use real physical units if it's applicable. Phase is measured either in radians or in degrees. These are VERY different things, so proper labeling with units helps to clarify what you are using. Similar if f is a frequency in Hz (or kHz or MHz) and T is a time (in nanoseconds or fortnights), it should be stated.
You could consider changing the Y-grid so that the grid lines match with meaningful phase values such as 2*pi (or an integer multiple thereof).
Along the same lines, you could alternatively normalize the phase to a meaningful value (like 2*pi) and keep the integer grid.
If the result was obtained by Fourier transform of a real valued time domain signal and the context is clear, consider omitting the left half (f<0) of the plot. The convention is to only plot positive frequencies and the complex conjugate symmetry is implied.

This may sound nit-picky, but making a good graph is typically time well spend. It really helps to clarify what's going on, simplifies communication and will help you to remember what you did when you look at it in 3 months from now.


Noise variance calculation for 16APSK Demodulator

This is to take into account the oversampling operation.

Symbol time = sample time implies no oversampling.

See AWGN model for more details about the conversion among EbN0, EsN0 and SNR.

  For example, if a complex baseband signal is oversampled by a factor of 4, then EsNo exceeds the corresponding SNR by $10\log_{10}(4)$.


How to transform flat-field corrected spectra to those under Standard illuminant E (ideal equal energy illuminant)

Is the flat field correction the only correction applied to the hyperspectral image? If the camera is not calibrated to radiance values then the fact that you need to apply shift/scale adjustments is normal and expected. The CIE XYZ color of a spectrum of light is defined in terms of radiance, the raw output of a uncalibrated camera is not radiance.

If the camera is radiance calibrated then here are a couple other ideas:

When you compare your computed color vs the known value are you scaling them both so that the CIE Y value of the white patch is 1 (or 100) ?

They need to both be exposure scaled equally to be compared, differences in the absolute scale will affect the CIEDE2000 values.

Also, you say the test target published colors are for illuminant E but was the hyperspectral image acquired under an approximation of illuminant E? If not the colors are not comparable.


How to plot magnitude response using $\tt freqz$ in MATLAB if I have a transfer function without a denominator

You simply specify a = 1 (because the denominator equals $1$). So you get

b = [1,sqrt(2),1];
a = 1;
N = 512;
[H,w] = freqz(b,a,N);

You can compare this to the analytical solution:

H2 = 1 + sqrt(2)*exp(-1i*w) + exp(-1i*2*w);
max(abs(H2-H))    % 8.0825e-16


The system is described by a linear difference equation with constant coefficients and as such, it is described in the same way as a linear time-invariant system. It is just the non-zero initial condition that makes the system non-linear, at least according to the common definition of linearity in system theory (homogeneity and additivity).

Such a system with non-zero initial conditions is also referred to as incrementally linear. An incrementally linear system responds linearly to changes in the input. For this reason it is very common, and it also makes a lot of sense, to discuss it in the same chapter as general discrete-time LTI systems.

As a final note, in practice we're mainly interested in stable systems. For such systems, the influence of non-zero initial conditions becomes negligible after some time, and, consequently, the system will practically behave as an LTI system after the transient caused by the non-zero initial conditions has died out.


How to transform large set of audio samples using FFT?

Can I transform the whole signal at once ?

Yes, this would be a 2400000 points FFT - something that is easy for a modern desktop computer (note that the sample rate is measured in Hz or samples/sec, not in Hz/s).

  If I cannot, is there away to combine the spectrum of different N sized segments?

Averaging the spectrum obtained on successive short windows is a common spectral estimation technique - but I'm not sure this is really what you need here.

  Is it reasonable to study the spectrum of a whole large signal at once ?

It makes sense if the signal is stationary. I can think of some applications in which this would be the case (for example you record 5 min of continuous hum/noise generated by a machine), but this would be an exception more than a rule - and I can't think of a practical reason why one would have to perform a FFT on such a large chunk of signal.

For signals whose characteristics change with time - such as speech or music - spectral analysis is more commonly performed on short segments. It all depends on what you want to do with the result of the FFT. If it is to extract a property of the signal (such as its pitch, or a metric like spectral centroid), make sure it is reasonable to assume that the property you want to extract is constant over the duration of the signal. If not, break down the signal into short segments during which this property can be considered constant.


How to find a know signal in another signal which contains other signals along with the signal i need in python

I made a little research on detecting the sync_start location on my received signal.

First I modulated the bits (sync_start bits + data bits + sync_end bits) with a sine signal as my carrier, with a frequency of $20 \text{ kHz}$ and a total length of 10000 bits (including both sync_start and sync_end).

I then converted this modulated wave into a .wav file
and played it on my PC, while recording simultaniously on an another PC, using Pyaudio (Python). That gave me a recorded .wav - my received signal.

Now as I knew my sync bits value, I modulated them on the same carrier wave and cross correlated this with the received wave.

Case 1:

$\text{bit length} = 10000$ (including sync_start and sync_end),

$\text{number of sync_start bits} = \text{sync_end} = 2287$

The plot of cross correlation is as shown :

Case 2:

$\text{bit length} = 15000$

case 3:

$\text{bit length} = 20000$

As the data length increases(keeping our sync_start and sync_end same in all cases), the exact time stamps of sync_start_mid index and sync_end_mid index gets more accurate.

Note : For data bits here, I took all of them as $1$'s -- that might be the reason for the smaller amplitudes in the middle. I will add random data and see if I am getting the same inference.

Please correct me if there is any mistake in my inference. As far as now I am now able to exactly pinpoint the synchronization between transmitted and received waves.

I also wanted to know the limitations of length of the sync bits (because when my sync bits length gets smaller, the dominant peak location cannot be found from cross correlation).

And I am using a BPSK scheme.


The meaning of correlating "Time Aligned" pseudo-random or walsh code sequences

Cross-correlation $r[l]$ between two sequences $x[n]$ and $y[n]$ is defined as:

$r[l] = \sum_n x[n]y[n+l]$.

You can say that $x[n]$ and $y[n]$ are "time-aligned" at lag $l^\star = \arg \max_l r[l]$. It is possible that the time-alignment is achieved at $l^\star=0$, which is the case you described above. However, in a communication channel with multipath delays, the alignment may happen at a lag $l^\star \neq 0$. Since you referenced DSSS receivers, I'd encourage you to read about a rake receiver and how signals from different fingers are combined.


detect to rising, stable and falling point in non-smooth rectangular wave

The usual approach to change detection is the CUSUM algorithm.

I've done an implementation that just addresses the level (mean) change issue. It's included (in R) below.

The black line is the noise-free data, the red line is the noisy data and the blue bars are the detected breaks (for this realization).

This just addresses the level change; to address the three zones: low level, changing level, and high level, you'll need to figure out the way to estimate the mean in the changing level (perhaps assume a rise time and fix the two levels?).

The code below is based on this document.

A slight change might get you the first and second changes easily: instead of working with the data, work with the difference of the data.  If I do that, then I get:

but I've had to add:

sigma <- 0.001
data_before_diff <- noiseless_data + rnorm(length(noiseless_data),0,sigma)
data <- diff(data_before_diff)
mu2diff <- 0.04
thresh <- 1

and then later:

 s[k] = (mu2diff - mu1)/sigma*(data[k] - (mu2diff + mu1)/2)

and

    tmp <- mu1
    mu1 <- mu2diff
    mu2diff <- tmp

To find the changes at the end of the signal, you'll probably need to change mu2diff to -mu2diff and check for that versus 0.

Note I've had to greatly reduce the noise variance and thus the threshold to get this to work.

R Code Below

# 30039
# CUSUM

N <- 30

mu1 <- 0
mu2 <- 1
noiseless_data <- c(rep(mu1,N), seq(mu1,mu2,1/N), rep(mu2,N*4), seq(mu2,mu1,-1/N), rep(mu1,N))

sigma <- 0.1
data <- noiseless_data + rnorm(length(noiseless_data),0,sigma)

thresh <- 100

breaks <- rep(0,N)
num_breaks <- 0

s <- rep(0,length(data))
capS <- rep(0,length(data))
G <- rep(0,length(data))
for (k in 1:length(data))
{
  s[k] = (mu2 - mu1)/sigma*(data[k] - (mu2+mu1)/2)
  if (k==1)
  {
    capS[k] = s[k]
  }
  else
  {
    capS[k] = capS[k-1] + s[k]
  }

  G[k] <- max(0,capS[k] - min(capS[1:k]))

  if (abs(G[k]) > thresh)
  {
    #capS[k] <- 0
    num_breaks <- num_breaks + 1
    breaks[num_breaks] <- which.min(capS[1:k])
    tmp <- mu1
    mu1 <- mu2
    mu2 <- tmp
  }
}

first_break <- min(breaks[1:num_breaks])
last_break <- max(breaks[1:num_breaks])

plot(data,col="red", type="l")
lines(noiseless_data)
lines(c(first_break,first_break),c(0,1),col="blue", lwd=10)
lines(c(last_break,last_break),c(0,1),col="blue",lwd=10)


How to separate the upwards-propagating from the downwards-propagating waves?

Instead of nulling the 2 quadrants, multiply them by $i$. Then the real part will be one image and the imaginary part the other image.

X = fftshift(fft2( x )); % x is the input (image above)
X( 1:floor(size(X,1)/2),    1:floor(size(X,2)/2)   ) = ...
    X( 1:floor(size(X,1)/2),    1:floor(size(X,2)/2)   ) * 1i;
X( ceil(size(X,1)/2)+1:end, ceil(size(X,2)/2)+1:end) = ...
    X( ceil(size(X,1)/2)+1:end, ceil(size(X,2)/2)+1:end) * 1i;
y = ifft2( ifftshift(X) );
imagesc(real(y)); pause;
imagesc(imag(y)); pause;

Since this is unitary operation the energy of this complex image will be the same as the original, as will the noise energy. Doesn't look very white to me though, seems there is differences between entire rows.

If you want to get rid of the artefacts, you could use a smoother function such as cosine of the angle of the spectrum:

[rows,cols] = size(x);
IM = fft2(x);
% Create polar co-ordinate of FFT
[ux, uy] = meshgrid(([1:cols]-(fix(cols/2)+1))/(cols-mod(cols,2)), ...
    ([1:rows]-(fix(rows/2)+1))/(rows-mod(rows,2)));
ux = ifftshift(ux);   % Quadrant shift to put 0 frequency at the corners
uy = ifftshift(uy);
th = atan2(uy,ux);
r = sqrt(ux.^2 + uy.^2);
% cos^2
R2 = cos(th +pi/4).^2 + 1i*cos(th +3*pi/4).^2;
R2(1,1) = 0;
I_a = real(ifft2(R2.*IM));
I_b = imag(ifft2(R2.*IM));
imagesc(I_a); pause;
imagesc(I_b); pause;

I don't think the energy remains constant in this case, but assuming white noise energy the noise of the resulting two images should be similar.

Using the polar form of the spectrum it will be easier to experiment and create smoother cutoffs. For example, you could 'narrow' the angle of the filter to isolate the wave further.


To compute the PSD, you'd need full amplitude/phase measurements of the noise signal.

If by "255 ac measurements", you mean measurements of the RMS amplitude by a digital multimeter, then you do not have the necessary phase information.

If by "255 ac measurements", you mean you are using an oscilloscope that is AC coupled, then you can consider the measurements to be a measurement of the voltage with the average value ("DC component") removed.   If we consider the voltage signal to be
$$
v_\text{total}(t) = v_\text{DC} + v_\text{AC}(t),
$$
$v_\text{DC}$ is going to be nearly entirely the "normal" contribution of the battery, and constant with time, and $v_\text{AC}(t)$ is going to vary with time and due entirely to thermal noise.  If your measurement is measuring $v_\text{AC}(t)$, then that is what you want anyways.

Once you have waveform samples, the PSD can be estimated using the periodogram (easiest, but noisy) or a more advanced method like Welch's or Blackman-Tukey method.


The best depends exclusively on the type of signal you want to modify, if you just need apply time scale in monophonic (Speech), PSOLA and TDHS can give to you very nice results, but It is not as easy as you read in books, find Pitch Marks (Glottal Closure Instant) in PSOLA is really a big problem, in TDHS if your pitch track are consistent you can scale the signal to the hell, for Time scaling unsig FFT (commonly known as phase vocoder) this works better in polyphonic signals and to get best results here do you need use a locked phase vocoder.

Find the length is a easy step, for stretch factor between 0.5 and 2.0 you can use some time domain technics: ola, sola, wsola, psola, tdhs, you can try keith lent (I consider a bit different of PSOLA) in the Lent code you can try 50% of crossfade hann window that works nice for my ears.


Use polyphase representation to construct a perfect reconstruction filter bank

This is a problem from Multirate Systems and Filter Banks by P. P. Vaidyanathan. In Section 4.6.2.E, there is a discussion on what he calls Euclidean Complementary Functions. Specifically, he discusses one of the results of Euclid's theorem:

  ... if $H_0(z)$ and $H_1(z)$ are relatively prime, there exists polynomials $F_0(z)$ and $F_1(z)$ such that
  $$H_0(z)F_0(z) + H_1(z)F_1(z) = c, c \neq 0.$$

In the text, he uses $H_k(z)$ and $F_k(z)$ to denote the analysis and synthesis filters, respectively. He uses $E_k(z)$ and $R_k(z)$ to denote the Type I and Type II polyphase components of the analysis and synthesis filter banks, so I'll do the same here.

If you take the polyphase decomposition of $H_0(z)$, we find that $E_0(z) = 1 + 3z^{-1} + z^{-2}$ and $E_1(z) = 2 + 2z^{-1}$. Because $E_0(z)$ and $E_1(z)$ are relatively prime, we can find an $R_0(z)$ and an $R_1(z)$ to fulfill

$$E_0(z)R_0(z) + E_1(z)R_1(z) = c, c \neq 0.$$
We can recognize that $c$ is the inevitable delay for $\hat{x}(n)$, so we should probably make $c = 1$. To cancel out the $z^{-1}$ terms so that $c = 1$, the order of $R_0(z)$ needs to be one less than that of $R_1(z)$, by inspection. So, we can represent $R_0(z) = \alpha$ and $R_1(z) = \beta + \gamma z^{-1}$, then solve for the coefficients. If we multiply everything out, we'll see that we want

$$\alpha + 2\beta = 1,\\
3\alpha + 2\beta + 2\gamma = 0, \text{and}\\
\alpha + 2\gamma = 0,$$
which is just solving

$$\begin{bmatrix}
1 & 2 & 0 \\
3 & 2 & 2 \\
1 & 0 & 2
\end{bmatrix}
\begin{bmatrix}
\alpha \\
\beta \\
\gamma
\end{bmatrix}
=
\begin{bmatrix}
1 \\
0 \\
0
\end{bmatrix}
$$
to get $\alpha = -1$, $\beta = 1$, and $\gamma = 0.5$. Now, we can solve the problem by composing $F_k(z)$ from the Type II polyphase components: $F_0(z) = z^{-1} R_0(z^2) + R_1(z^2)$ and $F_1(z) = z^{-1} R_0(z^2) - R_1(z^2)$.

There is actually also a nonzero scaling factor of 2 that comes from the DFT filter bank representation, because $\mathbf{W}^*\mathbf{W} = M\mathbf{I}$, where $M = 2$ because we have two channels.


Grayscale conversion before applying sobel gradient operator

The values used by OpenCV are values used by Luma coding in video systems. Your original code just used an equal weight for each channel.

This may or may not be important for your application (probably not really).


What's the state-of-the-art algorithms in image edge detection?

Despite its age, Canny Edge Detection is still a state of the art filter. The results produced by this algorithm make for it always being included in image editing software. Solid and descriptive edges that are often overly represented by other filters. It lacks the simplicity of, say, the Boolean Edge Detection, included in the paper, "Edge Detection Techniques: Evaluations and Comparisons" by Ehsan Nadernejad, Sara Sharifzadeh, Hamid Hassanpour.

The quality of Edge Detection and Segmentation filters are difficult to appraise, and algorithm speed is usually the most important consideration when implementing them.

For code and discussion, you should take your research to codeproject.com, and lookup Image Edge Detection.


Your problem is quite different than the vein extraction one. In fact, it is much mucg easier to solve with the "geometric" definition of the Hough Transform.

The Hough Transform is used exactly to detect lines in an image. It achieves this by integrating the brightness values of an image towards some angle $\theta$. So, for example, for angle $\theta = 0 ^ \circ$, the Hough Transform produces the sum of each row in the image, for angle $theta = 90 ^ \circ$, the Hough transform produces the sum of each column in the image...and so on.

The combined result of this is that points map to arcs, lines map to points and polygons map to specific configurations of points in the Hough Space.

The way to pick up a line in the image, in the Hough Space, is ultra simple. You just pick a maximum. But, that would give you an infinite line, not a line segment, which is what you have in your application.

The most straightforward way to pick up a line segment by the Hough Transform is to "mark" the pixels that are "responsible" for a particular accumulation. So, later on, when you pick up a maximum peak, you can "look" at the pixels that formed that sum and then connect them with a simple distance based rule. If distance is smaller than 4 pixels between any two pixels, then consider them as part of one line, otherwise create a new line.

For more information about all this, please see this link.

We now come to your case, you have two problems:

Detect line segments from GPS points
Connect line segments to form roads.

To solve #1, you can use the "geometric" definition of the Hough Transform WITHOUT having to work on the pixels of an image. Here is how this works:

Find the centroid of all GPS points. Call this the centre of your region.
Find the bounding box of all of your GPS points. These two steps define your "image space".
Run a Hough Transform with the center and over the bounding box of the points. As you are shooting "rays" down each different bearing, make sure that you also store the points that are responsible for the sum. The sum is simple: Shoot a ray down a direction, whenever it hits a point increase a counter and store its location.
Do a histogram of your Hough Space. It will help you in determining the threshold you have to apply so that you select the peaks that result in lines.
Apply the Threshold by selecting peaks that are above the thrshold value. You now have a bunch of infinite strong lines in your hands.
Look at the lists of points accumulated in each one of the directions and "integrate" them in lines if their distance is smaller than some threshold. You now have a bunch of line segments in your hands. These are your roads OR your road segments. The good thing about them is that they share a common point when they are connected (the junction). The bad thing about them is that they are unconnected when there is a "turn" or "curve".

You now have to solve problem number 2. Unfortunately, it is impossible to solve problem #2 without yet another application of a threshold. The philosophy here is simple:

Pick up a line, call it $a$.
Pick up all of the lines $a$ shares a common point with, call them $C$.
Calculate all angles between $a$ and the lines in the set $C$.
Merge two lines into one "road" if the angle is greater than some threshold. Obviously, if two lines form a $180^\circ$ angle, they are the same "road" but if they form a $90^\circ$ angle then they probably are not EXCEPT if set $C$ only has 1 element (1 line).

At the end of this process, you will have all connected segments that form roads.

It is not a difficult problem but it sure is a messy one.

Hope this helps.


Using Parks-Mcclellan (or others) for FIR differentiator?

There is no fundamental difference between the design of standard frequency selective filters (such as low pass, high pass, etc.) and differentiators. The only difference is that for a (linear phase) differentiator, the desired frequency response must be formulated as a purely imaginary function of frequency, unlike standard frequency selective filters, for which the desired response is purely real-valued. But that's a minor issue, and all important algorithmic details remain the same. Generally, algorithms only become different as soon as the desired response must be defined by a complex-valued function (with non-zero real and imaginary parts).

E.g., the Matlab function firpm, which implements the Parks McClellan method, can design FIR differentiators if you supply differentiator as filter type. The same is true for the least squares filter design function firls.

In sum, designing linear phase FIR differentiators is no more difficult than designing any other FIR filter.

Concerning low pass differentiators, take a look at this answer.


A pluck might produce significantly more broadband noise than a free string.  An FFT of such noise would show more relative energy outside of all the FFT result bins that are related to (F0 or overtones of) a single pitch.

Also, a free string has a more predictable decay rate in any FFT magnitudes related to the pitch across successive offset FFT frames windowing the sound data.  A sudden change from this decay rate (a stop or a big increase) might infer a disturbance, such as a pluck.


What audio effects, filters, distortions, etc. create a "vinyl" effect?

I don't know the precise formula for 1920's gramophone, but this is the general process for audio antiquing and should get you started. Depending on the exact settings, you can get vinyl, telephone, etc.

Convert to mono.
Add white noise.
Bandpass. I would start with a HP at 80 HZ and a steep LP at 8 kHz, but that's a guess. You'll have to use your ears until it sounds like the gramophone. You should try Chebychev and Butterworth filters if available.

You may also want to add "clicks" and "Pops". You can model those in a variety of ways. One way to start is with a Dirac function and filter it with bandpass functions. Depending on how strict your bandpass in step three is, you might be able to just add the clicks and pops as Diracs at the same time you add white noise, but then all your clicks and pops will sound about the same, which isn't very realistic.

Some other things you may want to try for more realism:

applying compression and saturation after step 1.
applying wow and flutter emulation either before step 2 or after (or, for most accuracy, two stages of white noise before and after the wow and flutter emulation).

You may want to ask on https://video.stackexchange.com/ to see exactly what filter settings they would recommend for 1920's gramophone in particular.

For tips on implementing basic audio filters, I have a blog post: http://blog.bjornroche.com/2012/08/basic-audio-eqs.html


How are two signals multiplied? And how is it different from convolving two signals?

Yes, you are correct. Multiplication in time domain means convolution in frequency domain and vice versa.  Multiplying your signals $x[n]$ and $y[n]$ will give an output:

\begin{align}
z[n]&=\{2\cdot 5, 4\cdot 1, 1\cdot 8\}\\ &= \{10, 4, 8\}\end{align}

Remember that this output is in time domain. When you convolve $x[n]$ and $y[n]$, you will get $z[n]$ in time domain as:

\begin{align}
z[n]&=\{5\cdot 2, 5\cdot 4+1\cdot 2, 5\cdot 1+1\cdot 4+8\cdot 2, 1\cdot 1+8\cdot 4, 8\cdot 1\}\\ &= \{10,22,25,33,8\}
\end{align}

Just flip one of the signals around zero and start moving right one place at a time. Multiply the corresponding points as you go along. The output has a larger sequence because convolution output has ${\rm length}(x)+{\rm length}(y)-1$ points.

To answer your query about where to use multiplication & convolution, assume you want to pass signal $x(n)$ through filter $y(n)$. The output of the filter $z(n)$ will be convolution of $x(n)$ and $y(n)$.
Now assume that you first converted from time to frequency domain, i.e. $X(e^{i\omega})$ and $Y(e^{i\omega})$ are frequency domain representation of $x(n)$ and $y(n)$. Now, to find $Z(e^{i\omega})$(output in frequency domain), you have to multiply $X(e^{i\omega})$ and $Y(e^{i\omega})$. To get the output in time domain i.e. $z(n)$ you have to apply inverse transform.

So that's how you use convolution and multiplication.


In which area the image dithering technique is used?

In image compression, when the number of bits per pixel are reduced, you might see an artifact called contouring.  Dithering is used to reduce this artifact.

For more information, please see

http://en.wikipedia.org/wiki/Dither


I don't know, but my guess on this is: just like air flowing through a flute, water flowing through a tube causes vibrations, tones. The strong tone with several overtones is pretty consistent with this being similar to an air instrument.


How multi layer map to different beams in LTE Beamforming

The answer is simply orthogonalization of radio resources. For any given OFDM symbol, exactly one user will be allocated certain sub carriers. This allocation is dynamic and informed to the user in what is known as the DCI, in the downlink control channel. Each user will decode data only on the subcarriers on which it has a grant allocated by the base station.

So each user is assigned a unique set of subcarriers for each OFDM symbol, hence when multiple users are multiplexed you do not get interference from other users. Because you are allocated a "unique orthogonal resource". Now if the base station wants to send multiple streams to the same user, then this is done by precoding matrices. 3GPP defines unique precoding matrices, the user estimates the channel and tells the base station about a probable precoding matrix and it is upto the base station to make the final call on the precoding matrix and inform it to the user.

Note: to your notion of beamforming on this question,

Beam forming, here is not in the sense of steering a beam in the direction of the user. But, precoding to transmit the stream in a given direction in the "vector space sense", where the vector space is formed by the N dimensional receiver space formed by the receive antennas. It is not a sense of physical direction as in conventional beamforming, it is more of beamforming in the receiver signal/vector space.


The only difference between the two transform results is a gain factor.  The transform should have the same energy as the original data sequence, so we can use that to check for the right result.  The energy of $X$ is-

$$
\sum X^2 = 484
$$
$$
\sum (5,9,11,3,−1,−1,2,0)^2 = 242
$$
$$
\sum (\sqrt{2}(5,9,11,3,−1,−1,2,0))^2 = 2\sum (5,9,11,3,−1,−1,2,0)^2 = 484
$$
This shows that of the two possible transforms, $\sqrt{2}(5,9,11,3,−1,−1,2,0)$ has the right energy and is thus the right choice.


Pwelch does not just use one fft, it applies several ffts , windowing, scaling and averaging.If you do all that you may get equality.  Here is example equality of fft & pwelch:

N = 94144;
x = randn(1,N);
pwr1 = pwelch(x,ones(1,N),0,N,1,'twosided');
y = fft(x);
pwr2 = abs(y).^2;
plot(pwr1);
hold
plot(pwr2*max(pwr1)/max(pwr2),'r--');


Spectral analysis of a cross correlation function

cross-correlation would determine the time-delay. but the scale of signals' length and the delay value should be considered precisely. I mean if you have two long signal of, for say, 100sec and you are looking for a time delay in an order of 1msec, you don't need to do cross-correlation over the whole 100sec. Instead, choose a proper window time (in msec order) of two signal and then perform the cross-correlation.


Reason for it might be the continuous adaptation, even in non speech regions. So, you can try Voice Activity Detector (VAD) for knowing the exact regions for adaptation of filter, and processing will be frame wise (20 ms).


Finding Accurate BPM and Beat Time Values in Audio

Comparing your results against data annotated by a human listener is the way to go. You can try downloading some of the evaluation data from MIREX or look for other research databases.


How to Check Separability of 2D Filter / Signal / Matrix

Nilesh Padhi, Welcome to the DSP Community.

The classic definition of separable means the data (2D) given by $ X \in \mathbb{R}^{m \times n} $ can be written as:

$$ X = \sigma u {v}^{T} $$

Where $ \sigma \in \mathbb{R} $, $ u \in \mathbb{R}^{m} $ and $ v \in \mathbb{R}^{n} $.
This is called Rank 1 Matrix.

How can you get those parameters and vectors given $ X $?
Well, the Singular Value Decomposition (SVD) is here to save the data.

The SVD of $ X $ is given by:

$$ X = U \Sigma {V}^{T} = \sum {\sigma}_{i} {u}_{i} {v}_{i}^{T} $$

You can see those match when $ {\sigma}_{j} = 0 $ for $ j \geq 2 $.
So what you should do is the following:

epsThr = 1e-7;
[mU, mD, mV] = svd(mX);
vD = diag(mD);
if(all(vD(2:end) < epsThr))
    vU = mU(:, 1);
    vV = mV(:, 1);
end

We checked if indeed the singular value of 2 and onward are small.
If they do (You can decide to what degree of small by epsThr) then it is separable and the vectors are vU and vV.

In your case:

mX = [1, 3; 2, 6];
[mU, mD, mV] = svd(mX);

vD = diag(mD);
disp(vD);

The result is:

vD =

    7.0711
    0.0000

Since vD values are zero besides the first element (Only single non vanishing Singular Value) it is separable.

Indeed you can see that:

mD(1) * mU(:, 1) * mV(:, 1).'

ans =

    1.0000    3.0000
    2.0000    6.0000

As expected.

This method is really useful in Image Processing when we want to convolve with 2D kernel and we find it is separable and hence we can apply the 2D convolution using 2 1D convolutions (Along Columns / Rows).

In that case we define $ \hat{u} = \sqrt{{\sigma}_{1}} u $ and $ \hat{v} = \sqrt{{\sigma}_{1}} v $ where $ {\sigma}_{1} $ is the Singular Value.
Then we convolve $ \hat{u} $ along columns and $ \hat{v}^{T} $ along rows.


FFT function isn't working on imported signals greater than a second

Here a short example of how to deal with FFT in octave/matlab:

function freq(npts, len, basefreq)
    x = [0 : npts - 1] / npts * len;
    y = sin(basefreq * x * 2 * pi);
    F = abs(fft(y));
    half = npts / 2;
    F = F(1 : half);
    plot([0 : half - 1] / len, F);
    printf("Base frequency = %g\n", (find(F == max(F)) - 1) / len);
endfunction

try this function with different parameters.

Here npts is size of data vector; len - time length of data (in seconds); basefreq - base frequency for sine function.
First we fit a time axe (x), after that we calculate sinusoid with given frequency, get modulus of its FFT, trim it to half (BTW, another variant is not to trim, but to move frequency zero to middle of fftshift'ed data). After all we plot result of FFT and print out founded base frequency.


Image Brightness on 2D Fourier Transformation

If I understand correctly only phase would change not frequency. Is that correct?

Assuming "increasing brightness by 10%" just means multiplying each pixel's value with 1.1:

No, the phase won't change either. Since the Fourier Transform is a linear operation, you just multiply every point in the spectrum by the same factor. That doesn't affect phase, since 1.1 is a real number.


Why do we quantize high frequencies coarsely, shouldn't it be the opposite?

In Huffman coding frequency refers to the number of times a particular sequence of discrete source symbols has appeared in the input. It is not the same concept as frequency in the continuous space such as sound frequency. As such, it has nothing to do with psychoacoustics or the cochlea.

We use fewer bits for the higher frequency source symbols because frequent source symbols have lower entropy and assigning fewer bits to low entropy symbols gives better compression.


I presume you are correct in that the small difference in the original signal fft's is due to leakage. However, the huge difference in the windowed results arises because you simply do not have enough data. Your original signals only have data for one period of the sinusoid. When you window this, you are clamping the edges of your signal to zero, effectively leaving you with only the middle portion of the signal. So you are left with signals that have only a fraction of a single period worth of information, definitely not enough to have any accuracy in the extracted frequency analysis. To illustrate this, here are matlab images of your original functions before and after windowing:

and you can easily see the loss of much of the meaningful frequency information from the signals. A better candidate for windowing would be something like below, using 10 periods of your sinusoid signals:

which demonstrates the balance between retaining frequency information and clamping the edges to reduce leakage. This paper provides a more in depth look at different windowing functions and their performance on different classes of signals if you are interested.


Problems involving the derivations of CTFT $X(j\omega)$ for a given signal $x(t)$ or the reverse: finding $x(t)$ for a given $X(j\omega)$ can be solved in a number of ways. Fundamentally there exists the so called fourier integrals: the Analysis $$X(j\omega) = \int_{-\infty}^{\infty} {x(t) e^{-j\omega t} dt }$$
and the synthesis: $$x(t) = \frac {1}{2\pi} \int_{-\infty}^{\infty} {X(j\omega) e^{j\omega t} d\omega }$$
to directly compute what is being asked. But due to several reasons we also develop practical, shortway techniques to find those signals, as a very natural extension of engineering practice.

One such practical method of finding the forward or inverse fourier transforms involves the systematic application of known pairs with partial fraction expansions. Which enables the bypass of the integrals and yields fast answers.

First of all you shall apply the method of partial fraction expansion to your given CTFT $$X(j\omega) = \frac {3}{(1+j\omega)(2-j\omega)} = \frac {1}{1+j\omega} + \frac {1}{2-j\omega} $$

Then you shall consult into a table of fourier pairs, which would guide you to find which signal has transforms of this kind. From those standard tables one can see that $x(t) =  e^{-at} u(t) $ has the transform of $1/(a + j\omega)$ and also $x(t) =  e^{at} u(-t) $ has the transform of $1/(a - j\omega)$, for $ a> 0$

Henceforth, you would conlcude by finally invoking the linearity property of CTFT that $$ x(t) = e^{-t} u(t) + e^{2t} u(-t)$$


How to determine relative distance between 3 in-phase transmitters by analysis of received signal

Check this out, it does calculate the similarities of 2 or more signals
http://www.mathworks.com/help/signal/examples/measuring-signal-similarities.html


Digital PLL loses lock every time mesage crosses zero

You are having trouble because that's not a Costas loop.  A Costas loop uses demodulated data in some form to change the phase that's expected from the signal.

You're just taking the I/Q demodulated signal and applying it to the atan2 function; that makes a sort of linearized extended phase detector, but without determining that the phase should have flipped by 180 degrees, your loop is naturally trying to follow the phase that's actually coming in.

You could make your loop into a simple Costas loop by changing the atan2 function into a multiply, and following it by a proper loop filter.


Software-Emulated Amplifier: Cracking while change Volume

I’ve heard it called popcorn noise before.  Think of the output signal as the product of the input signal and the volume signal. A discontinuity in the volume signal causes a discontinuity in the output signal. The solution is to limit the bandwidth of the volume signal.  The most common implementation I’ve seen is to limit how fast the volume can change.  Another option would be to feed the volume signal through a low pass filter.  Note that the output of an LFO is already (pretty much) band limited, which explains why you wouldn’t hear any noise in that case.


Median Filter for salt and pepper noise removal

Median filter is considered good because unlike averaging filter which ruins the edges of an image by blurring it to remove the noise, median filter removes only the noise without disturbing the edges.

Well, median filter is the best and only filter to remove salt and pepper noise.

Hope this helps:)

Thank you!!!


Confusion Related to the Laplacian Filter

The Laplace filter is the sum of the differences of neighbourhood pixels compared to the central pixel. Its output will be zero in regions where the intensity is constant, and non-zero when there is a transition in intensity (i.e. an edge).

For simplicity's sake, consider what happens in one axis only, where the filter kernel would be [1 -2 1]. If we have a sequence of values with a step transition then the input and output will look like this:

input:  10   10   10   10   10  100  100  100  100  100

output:       0    0    0   90  -90    0    0    0


How to remove/reduce noise by using the Morlet wavelet?

I'm not sure of the reliability of your source... I worked with wavelets on my master's thesis and the basic wavelet transform is a convolution of the wavelet an your signal.
In most practical cases, it's best to go through Fast Fourier Transform (FFT), as it's mathematically equivalent.

With WT being the wavelet transform, and FFT() being the Fast Fourier Transform of what's inside the brackets, iFFT() is the inverse-FFT. You would compute the WT as such:

FFT(WT) = FFT(signal) * FFT(wavelet)
WT = iFFT(WT)

Now, for the wavelet formula itself, I believe there has to be a specific scaling factor to get a reliable result. Here's what we used:

$$\Psi = exp\Big(-\frac{ik_0x}{\lambda} - \frac{x^2}{2\lambda^2}\Big)$$

Where $\lambda$ is the wavelength (or scale) and $k_0 = 5.336$ (source: J.F. Kirby, Computers & Geosciences 31 (2005) 846–864)
I did not include the normalization factor because you just have to divide your wavelet by it's highest value, and subtract its average, and you get the same result with less headhache.
Off course, you can interchange $1/\lambda$ by $2\pi*frequency$ in the formula.

In order to qualify for the "Morlet" name, I think you have to use the same wavelength in both places in the formula. Else, you over-sample or under-sample the frequency part of your signal. $k_0$ is the best ratio between the two, as explained by Kirby in his article.


I'm still not sure that I understand what your question is, but this specific one is easy;

  My question is why do I need the inverse of the integral of that window?

The issue is that you'll apply a window to the data $x$:
$$
x_2(t) = x(t) \cdot w(t)
$$
that means that the sum of $x_2$ is going to be different from the sum of the (rectangular windowed) $x$.  All the inverse of the integral is doing is rescaling things to attempt to make the sums equal.

As far as your other question goes:

  Are my units going to be nT?

That's all down to calibration. You have millvolts.  Your expression for $H(\omega)$ has a calibration factor $0.8V/nT$. To me, that means that 0.8V (or 800 mV) is 1 nT.

Tesla are one kilogram per second squared per ampere in base SI units.  Usually, with any measurements calibration is what equates the voltage (or current, usually) measured to the units of you're actually interested in.

All $H(\omega)$ seems to be doing is applying a bandpass filter to your data.


Signal Estimation after detection Part 2

After signal detection, how to estimate the clean signal $s(t)$?

Matched filtering is used to detect the presence of a known signal in noise. There is no estimation part when you are talking about a matched filter. The estimate part comes after you have done the matched filter and need to estimate the symbols.
It looks like you are talking about a communication system context, but the matched filtering in those systems is not done by matched filtering the symbols as in your example code. For each symbol, the true symbol is unknown and the pulse shape is known, and that is why the matched filter uses the pulse shape to matched filter. Perhaps all that you want to know about matched filtering can be found here, https://dsp.stackexchange.com/a/9389/31316, and you should give it a read.
Edit
The OP questions about the zero threshold and whether there are other methods. That threshold is set through using the maximum likelihood solution, or minimum distance decoding. For the case of BPSK, you want to decide if $x$ contains $s_1=+1$ or $s_2=-1$, and the ML rule is to choose $\hat{s}=s_1$ if:
$$ |x-s_1|<|x-s_2| $$
This can also be interpreted as partitioning the IQ diagram into different regions (Voronoi regions), a region for $s_1$ and a different region for $s_2$. These regions end up looking like this (note how for BPSK the threshold line is at zero):

For full derivation, see Tse's Wireless Comms book section A.2.1, https://web.stanford.edu/~dntse/Chapters_PDF/Fundamentals_Wireless_Communication_AppendixA.pdf.


dB SPL readings differ from external sound level meter?

Your calculations are correct, although you should do some kind of sliding average over a few hundred samples or so. The problem seems to be that the sound level meter employs an A-weighting, while your code does not. If there is low frequency noise in your measurement environment, the readings at "silence" can easily differ by 40-50dB. Try measuring a test tone at 1kHz, results should be closer together. Also, have a look at

https://en.wikipedia.org/wiki/A-weighting

You can see the huge low frequency differences there.


Transmission capacity in "duobinary coding"

If you want to transmit a symbol sequence $A_k$ using baseband pulse amplitude modulation (PAM), the transmitted signal is

$$s(t)=\sum_kA_kp(t-kT)$$

where $p(t)$ is the transmit pulse, and $T$ is the symbol interval. If you want to avoid intersymbol interference (ISI), the pulse function $p(t)$ must satisfy the Nyquist criterion, which says that its value must be zero at multiples of the symbol interval:

$$p(0)=1,\quad p(kT)=0,\quad k\neq 0\tag{1}$$

The minimum bandwidth required to satisfy (1) is $1/2T$, but the corresponding pulse $p(t)$ is an ideal low-pass filter, i.e. a sinc function, which is of course impractical. So in practice, Nyquist pulses with a higher bandwidth (excess bandwidth) are used. This means that for a given bandwidth, the possible symbol rate is reduced as compared to the theoretically maximum possible value.

With duobinary signalling practical systems can be built with zero excess bandwidth, because the data symbols are filtered such that a spectral zero is introduced at half the baud rate. Note that this is a very simple filter unlike the ideal low-pass filter required for the minimum bandwidth Nyquist pulse. The price paid is the deliberate introduction of ISI and the introduction of a ternary (instead of binary) signal, which requires a higher SNR to achieve a given error probability.


A register exchange method implementation of Viterbi decoder avoids trace back and a nominal initial delay (one time delay or boot up delay) is unavoidable. But this method requires more memory resources.


Meaning of frequency and bandwidth of a signal, despite the fact that we do not know the signal

Some signals are generated by processes that are physically capable of generating frequencies only in a specific frequency range. For example, the human vocal tract can only generate signals between roughly 50 and 10,000 Hz. This does not mean that every human can achieve that range, or that there are people who can generate frequencies slightly above or below. When it is said that human voice has a certain bandwidth, it does not mean that every instance of human voice will actually cover that entire range; it only means that you will not often find frequencies outside that range.
Other signals are limited in bandwidth by design. One example is voice in analog telephony; the voice is filtered to a range of approximately 300 to 3400 Hz, and it is impossible to find frequencies outside that range in the system. Another example is analog television: the signal is designed to have a bandwidth of 6 MHz, and except for equipment malfunction, you'll never find a TV signal with frequencies outside that range.


Continuity and its relationship with asymptotic spectral decay

I suppose it would be sufficient to show that a discontinuous function in the time domain decays as $1/\omega$ in the frequency domain (e.g., $\textrm{rect}\Leftrightarrow\textrm{sinc}$), and then use the fact that each integration in the time domain adds a multiplicative factor of $1/j\omega$ in the frequency domain.