What is the underlying concept behind Bandpass sampling?

Yes, the signal is changing "every nanosecond" for any non DC signal.  But if the bandwidth is really limited to 0.001 of the carrier frequency, then those changes are very predictable, as the changes result in a carrier waveform that is a nearly perfect sinusoid of a known frequency.  The signal can't differ too much from that highly predictable sinusoidal waveform unless the bandwidth is increased.  Assuming the bandwidth doesn't increase, any changes to the near perfect sinusoid are so small that one only needs to check the waveform every few hundred cycles of the sinusoid to see if it has drifted a bit from the expected point on the expected cycle.  Thus under-sampling at somewhere above 0.002 cf can tell you all you need to know in addition to knowing that the bandwidth was absolutely constrained to 0.001 cf or below.

But the sampling jitter requirement has to be on the order of a very tiny fraction of a nanosecond, or else one will completely miss the expected point on the expected cycle of the 1 Ghz sinusoid.


scipy.signal.correlate2d(A,A) and scipy.signal.fftconvolve(A,A[::-1, ::-1]) will give the same result. The second one is much, much faster.


"So do I only use the sinc functions for the samples that are closest to my x value?"
Yes, when you are truncating.
"Is this what is meant by "windowed sinc"?"
Yes.  The sinc goes to infinity, calculating that is impractical.
"How many samples should I go in each direction away from my x value?"
The sinc diminishes the further away you are.  At some point additional points become insignificant (no wider than a whisker on a gnat).
"Additionally, my DAW has something called "32 point sinc" resampling. What is the "32 point" supposed to mean in this case?"
I'm pretty sure it means a 32 point window, centered on your current sample.
Sinc interpolation is one form of interpolation.  It is the "Fourier compatible" one as it gives limits on a bandwidth basis.
Here is another link that may be even more useful:  Multi-channel audio upsampling interpolation
On re-reading this reference, Olli has already done this analysis, spectacularly.


breaking the input signal into segments while calculating DFT

How they break up the input signal into 256 segments?

This is usually done by a sliding a window of a given width (here 256). Depending on the step-size this results in overlapping (if smaller than 256) or non-overlapping (if equal or larger than 256) data segments.

  How each segment was multiplied with the Hamming window?

Usually, by an element-wise multiplication.

  How they calculate 256 DFT points?

The number of DFT points is, if I am not mistaken, determined by the length of the input signal.

  How the frequency spectra was averaged?

All resulting spectra were summed up and divided by the number of spectra.


ICY feature detector and edge detector in opencv

I would suggest to smooth the image before you apply the edge detection.
Apply bilateral filter on the image before the edge detection( cv2.biltralfilter).


Squaring signal for amplitude demodulation

If $m(t)$ is the message signal with $|m(t)|<1$, the AM signal can be written as

$$s(t)=(1+m(t))\cos(\omega_ct)\tag{1}$$

where $\omega_c$ is the carrier frequency in radians per second. Squaring $s(t)$ gives

$$s^2(t)=(1+m(t))^2\cos^2(\omega_ct)=(1+m(t))^2\frac12(1+\cos(2\omega_ct))\tag{2}$$

Low pass filtering $(2)$ gives

$$x(t)=\frac12(1+m(t))^2\tag{3}$$

from which the message $m(t)$ can be retrieved:

$$m(t)=\sqrt{2x(t)}-1\tag{4}$$


How I can perform SIMO channel convolution in MATLAB

Your question is specific to how to do the convolution with a matrix multiplication in the case of a SIMO system.

Below is a graphic from Compensating Loudspeaker frequency response in an audio signal showing how the Toeplitz matrix is used to do SISO convolution with matrix multiplication:

In your case "r" in the graphic is "x" in your code, and the Toeplitz Matrix A is formed by placing x in the column and shifting this along the diagonal over the number of columns given by your channel length (in your case 12). So if your transmit signal was 160 samples long and your channel 16 long, A would be a 175 x 16 matrix. As shown in the figure above, if you only had a single channel, the channel would be multiplied as a column vector resulting in the matrix product of 175x16 with 16x1 resulting in a solution as a column vector of 175 x 1. (Remember to do the matrix multiplication, the inside dimensions must match and the result will be the outside dimensions).

For multiple channels as your case, simply multiply by the 16x4 channel matrix instead. This will result in 175x16 mulitplied with 16x4 resulting in a matrix result of 175 x 4, representing a column for each received signal.


Tips for programming scientific articles

Unfortunately, there is accumulating evidence on lack of reproducibality of results across different fields (for example, for another example, for yet another example).

As a response to that, professional bodies and research teams are publishing guidelines for reproducible research and you might have to search for something like that in your field. It would be great to cite something along the lines of "...in our work we have tried to take on board the Reproducible Research guidelines published by...[citation] wherever it was possible".

Some of those are as follows:

Ten simple rules for reproducible computational research
A guide to reproducible code in ecology and evolution
Recommendations to funding agencies for supporting reproducible research by the American Statistical Society
Reproducible Research for Scientific Computing: Tools and Strategies for Changing the Culture
The ROpenSci project's Reproducibility in Science. The specific reason I am mentioning this is because it has a very nice "Further Reading" section, albeit specifically addressing R users.

And of course, on the wake of these guidelines, there are projects like the Open Science Framework which try to provide implementations to guidelines and propagate best practices.

While you can read lots of articles, about these things, the point is to create such conditions that demonstrate clearly how you arrived to your results. And this has to be done in such a way that it does not require your input. You don't have to be there to explain anything (this is also what the scientific papers are supposed to be about).

This is very powerful for three reasons:

Obviously, if there is any doubt by anyone about your results, they can simply try to repeat everything themselves or, if possible, simply download a virtual machine with "the whole lab" already set up and go through your whole process.
If you have made some sort of mistake, it is easy for someone to point it out and it is also a way to demonstrate that this research was performed with the best knowledge at the time it was conducted.

If you have made some sort of mistake and someone else spots it everybody wins.

If you cannot reproduce the results then you simply say "Here is the whole setup, here are all of my workflows and data, here is how the experiments were setup (if you cannot share the data directly). I have tried my best to reproduce these results but so far I have not been able to do so".

This last point is probably more important for what you are trying to do because sometimes, there is a lot of frustration as to why "our numbers don't satisfactorily match XYZ's numbers".

Hope this helps.


Separating 'signal' from periodic 'noise'

If you just need to detect when a gaussian peak occurs, you could try this:

Use the scipy hilbert function to determine the amplitude of the signal.

Pass the result through a high pass to remove the DC component.

Pass the result through a low pass to remove high frequency noise.

Use a threshhold function to detect peaks.

This will tell you when a peak occurred, and should be good at rejecting constant frequency and amplitude signals, as well as rejecting white noise of more or less constant amplitude.

It won't necessarily give you the shape of the peak. The peak you will see in a plot of the result (before the threshholding) may well resemble the original peak, but it won't be a perfectly accurate representation of it.

The values of the high pass and low pass filters can be juggled to get rid of noise, and to provide an appropriate pass band for the signal you are looking for.


In response to MBaz's comment, if we define the impulse response (edit) for Linear Time Invariant systems as

  "the function $h(t)$ such that the mapping $x(t) \mapsto y(t)$ is computed as $y(t) = x(t) * h(t)$"

then notice that if $x(t) = \delta(t)$ you get $y(t) = \delta(t) * h(t) = h(t)$.

Thus louism's statement a) and MBaz's comment agree.


why is the signal divided into epochs in EEG classification?

Basically, to perform a spectral analysis. Take a look to this: http://support.ircam.fr/docs/AudioSculpt/3.0/co/Window%20Size.html
I do not know if you refer to some specific approach, but usually, the main EEG classification approaches use frequency analysis (also combined with temporal analysis) to extract the features that will be used to train the classifier or infer the prediction.
Cheers


Salt and Pepper impulse Denoising opencv

Can we consider these noises as salt and pepper noise.? Is there something else that I am missing?

several pixels getting erased to either zero or one -> yeah, that fits my definition of salt and pepper noise. However, there might be other definitions around, including assumptions of noise autocorrelation, and whether or not your noise fulfills these is impossible to answer without knowing what definition you want to follow.

Anyway, yeah, let's call this salt and pepper, doesn't hurt anyone, does it?

  What could be the better suggested method? Currently, I am planning to implement Switching Median filter by Pei-Eng et al, but I would like to know if this could be the right track.

I'd start by looking for pixels that are either exactly white or exactly black – those are your only candidates for noise pixels. Leave the rest alone.

Among these pixels, I'd look at the 8 neighbor pixels. If 3 or more of them are a similar intensity, keep the pixel. It's probably properly white or black.

For the other pixels, a good estimate is probably an average (instead of a median) of the 8 neighbors plus itself, weighing itself maybe a bit less. Generally, a continuation of the structure of the surroundings.

By the way, be careful not to be too smart – the observer of the denoised picture doesn't know you've denoised it, and wrong diagnoses based on overcorrected imagery are thoroughly possibly. (I've lost a bit of faith in the ability of different physicians to consistently diagnose based on the same MRI, but coming from a machine they don't know, after having done a signal processing project in medical diagnosis.)


Is it correct to subtract a low-pass filtered signal from the original signal and use the result as a "high-pass"?

In theory you can do this, but in practice it is difficult to do because the time and phase alignment must be pretty good for it to work.  If the alignment is good you will get the destructive interference that you are seeking.  If they aren't, you will get some constructive interference.  Even worse, whether they are destructively or constructively interfering will depend on the frequency- i.e. you can get both constructive and destructive interference at the same time.  It can work, though, if you are only filtering out fairly low frequencies, since their timing requirements are the loosest because they change so slowly.

Short story- it's possible to do but is difficult enough that it generally makes sense to just do a high-pass filter.

A relatively straightforward way to create a bandpass filter is to create a low-pass filter, and then modulate it to the center frequency you want by multiplying it with a sinusoid of that frequency.


Problems Using FFT to Compute Impedance in a Model Neuron

I don't know how much background you have with FFT, so I apologize if I list a few things you already know.

In general, I think people have a tendency to throw an FFT at a problem and hope it works, but you may need to consider what the transform is really all about. First, when you take an FFT, it's implicit in the calculation that your signal is periodic. So, even though you're feeding it a finite number of samples with frequency content only between f0 and ff, the calculation makes the assumption that before time zero and after your samples, the waveform repeats forever. If your waveform's beginning and end don't match perfectly, there's going to be a lot of frequencies you wouldn't expect in your FFT. Basically, it's like there's a discontinuity. So, often, you'll window the signal to force it to zero at both ends so each periodic clone of your waveform meets the next smoothly. This obviously also has consequences for the frequency content of the FFT, so you need to do this carefully also.

Second, if your waveforms have any DC content (a constant offset, frequency = 0 Hz), it usually makes sense to remove it before your FFT unless DC is of interest to you. Often, low frequencies are huge and swamp the stuff you care about. Doing some high-pass filtering with a very low cutoff will help you see the interesting content of your plot more easily.

Third, when you view your plot on a linear frequency axis, you'll see that there are Fourier coefficients for both positive and negative frequencies, and the plot is symmetrical about 0 Hz. For qualitative assessments, you probably only need to look at the positive half.

Fourth, consider that you're taking the ratio of output over input, and anywhere your input FFT is near zero, that's going to blow up, and may make it hard to see what's going on.

So, I'd say, do some preprocessing of your signal to remove things you don't care about, and to prevent edge effects. Then post another plot and let's see how it looks. Since you're simulating the input and output here, there may not be a need, but in practice, I'll often use Welch's method to compute power spectral densities, and then compute an H1 or H2 estimation of the frequency response. Here, we can hopefully get a result with a simple analysis.


A causal discrete-time LTI system is marginally stable if none of its poles has a radius greater than $1$, and if it has one or more distinct poles with radius $1$. So a system with poles at $z=1$ and $z=-1$ is marginally stable (if there are no other poles outside the unit circle).

A causal discrete-time system with all its poles strictly inside the unit circle is called asymptotically stable. If at least one of the poles is outside the unit circle, the system is unstable.

For anti-causal systems just replace 'outside' by 'inside' and 'greater than $1$' by 'smaller than $1$' in the above definitions.

For non-causal systems with two-sided impulse responses the region of convergence (ROC) is an annulus centered at $z=0$. Such a system is asymptotically stable if the ROC contains the unit circle. If the ROC is limited by the unit circle and if one or more distinct poles are on the unit circle, then the system is marginally stable.

In sum, we get the following ROCs for marginally stable systems:

causal: $1<|z|<\infty$
anti-causal: $0\le |z|<1$
non-causal (two-sided): $1<|z|<R,\; (R>1)\quad\text{or}\quad R<|z|<1,\; (R<1)$

where in all cases it is assumed that all poles on the unit circle are distinct.


Which Fourier Transform properties should be applied first when they are combined together?

If you do things right the order doesn't matter. Note that replacing $n$ by $n+4$ in $x[-n]$ results in $x[-(n+4)]=x[-n-4]$ and not in $x[-n+4]$.


Convolution:

$$ y(t) = h(t) \circledast x(t) = \int\limits_{-\infty}^{\infty} h(u) \, x(t-u) \ \mathrm{d}u $$

Cross Correlation:

$$ R_{xy}(t) = \int\limits_{-\infty}^{\infty} y(u) \, x(t+u) \ \mathrm{d}u $$

The difference between the two is effectively the sign on $u$ in $x(t-u)$ in the integral.  That correlation is like convolution but with one of the signals flipped left to right is essentially the basis of the Matched Filter.


How to Clean Black and White Images Which Are Results of Edge Detection

Turns out, the trick to solve this was just to realize that from frame to frame the background was changing, but the text was staying constant. So by doing an | operation on the pixel values of a few consecutive frames, it was crazy simple to get only the text.


The general formula for the exponential rise in frequency is given by

$f=A+B(1-\exp(-\alpha t))$

We need to solve for the values of $A$, $B$ and $\alpha$. Let's assume you want an initial frequency of $20$ at time $t=0$, and an asymptotic frequency of $200$ Hz at time $t=\infty$.

Substituting $t=0$ and $f=20$ into the above equation we find $A=20$. If we substitute $t=\infty$ and $f=200$ we find that $B=180$. For $\alpha$ you really have a free choice - it controls how fast the frequency rises. If you have another control point say at time $t=2$ we want the frequency to be at $95%$ of the final frequency i.e. $f=.95 \times 200=190$ then we can substitute these into the above equation (Since we know what $A$, and $B$ are).

$190=20+180(1-\exp(-2\alpha))$

Solving gives $\alpha= -0.5\ln(10/180)=-1.445$. If you want a faster rise then you'll need to choose a larger value of $\alpha$.

To find the expression for phase we simply integrate with respect to time to give

$\phi =2\pi(At+Bt+\frac{B}{\alpha}\exp(-\alpha t))$


Why isn't this AM Demod GNU Radio flow graph working with different sampling rates?

In order to make the filter calculations easer, you should specify its sampling rate equal to the source's. Assume that:

$D$ is the filter's decimation factor
$f_s$ is the source sampling rate
$f_c$ is the filter's specified cutoff frequency
$f_f$ is the filter's specified sampling rate

The actual cutoff frequency from the point of view of the filter's output, $f_a$, is calculated as: $$f_a=\frac{2Df_c}{f_f}\frac{f_s}{2D}=\frac{f_cf_s}{f_f}.$$ This is because the filter doesn't care what the actual sampling rate is; what it needs to know is the cutoff frequency as a fraction of the normalized Nyquist frequency, which is $2Df_c/f_f$ after accounting for the decimation.

In your first flowgraph, the actual cutoff frequency is $5\cdot256/64=20$ kHz, not 5 kHz (the same problem occurs in the second flowgraph; in the third, the cutoff is 25 kHz). If you specify the filter sampling rate equal to the source rate ($f_s/f_f=1$), then your signal will actually be filtered at 5 kHz.

You also need to make sure that the product of the source sampling rate, times all your interpolation, divided by all your decimation, is equal to the sink sampling rate. In the third flowgraph, the sampling rate into the sink is 16 kHz, not 48 kHz.

Having said that, I think there's another problem in your flowgraph. I assume that you're demodulating AM DSB-LC. If the RTL-SDR device is centered exactly at the carrier frequency, then you only need the low-pass filter; no complex to mag conversion is necessary. You only need to keep the real part of the signal (send the imaginary part to a null sink -- you can do this right after the RTL-SDR and save some computations).

You only need to do envelope detection when the AM signal is not in baseband. In that case, you calculate its absolute value and then low-pass filter it. You also need a final high-pass filter to eliminate the DC bias present at the LPF output.


What does it mean to train and test a data for feature extraction?

This is a concept in supervised machine learning.

Train data: Used to train your supervised ML model. This data contains both the input and the desired output, which is compared with the output from your model. Your model is hence trained to reduce the prediction error or cost function.

Test data: The data for which you want to predict the outcome. The desired output is unknown. Based on the training data, your model predicts outcome for the new test data.

You should become more familiar with the concepts of machine learning first, before asking questions here. There is a plethora of resources on machine learning.


Minimum sample frequency of IMU accelerometer and gyroscope

I'm not sure I can state anything better than Marcus's or Maxtron's comments (it's that way because it's always been that way / because the hardware is limited to that rate), but another way to think about it is:

The Nyquist rate is 62.5 (=125/2) Hz.
Can your acceleration or rotational angles (the things measured by the accelerometer and gyroscope) change at a frequency faster than that?

If your instrumented system cannot, then this sampling rate is OK. If your system can move faster, then the sampling will miss something.


Do I need an overlapping in STFT if I want to play the signal back?

Note that the algorithm in the article uses a windowed FFT for analysis.  Window functions help avoid certain FFT artifacts, but they strongly attenuate any signal at the ends of the window.  So if you don't overlap, almost all the audio (any pitch changes, etc.) at both edges of each window would be lost, and you would end up with a choppy re-synthesized result.  Also, the window function usually peaks in the middle of the window, and is thus produces more accurate analysis there.  So the closer together the window centers, the more accurate the analysis involved in the pitch shift algorithm.  But the window needs to be longer than this distance for accuracy in frequency estimation.  Thus the overlap of 75% instead of just 50% in the article.

Synthesis for playback often involves the phase vocoder algorithm plus cross-fading of overlapped windows to avoid discontinuity ticks and other glitches at the synthesis window boundaries.


Auto-covariance of the product of deterministic and wide-sense stationary signal

Let me give a clue
$$ E\{ y \} = E\{ x g \} = E\{ x \} g
$$


ADC performance simulation: How to calculate SINAD from FFT?

If SINAD can be determined from only expected value and variance then it is possible to determine how SINAD transforms. Variance $\sigma^2$ is preserved whereas expected value $\mu$ grows as $ \sqrt{N} $ where N is the sampling set size. Noise + distorsion is assumed to have variance $\sigma^2$.

Thus the SINAD value would then be determined to become
$$
{\mathrm {SINAD}}={\frac {P_{{\mathrm {signal}}}+P_{{\mathrm {noise}}}+P_{{\mathrm {distortion}}}}{P_{{\mathrm {noise}}}+P_{{\mathrm {distortion}}}}} =
{\frac {N |\mu|^2+\sigma^2}{\sigma^2}} = N \cdot SNR + 1
$$

I can explain this in more detail if needed.


Why is the median filter called 'edge-preserving'?

Non-linearity

A linear filter is mathematically described by the convolution sum (for discrete signals) and the convolution integral for continuous signals. The median cannot be found using a linear function except in the trivial case where you have a discrete filter of size 1, which is why the median filter is non-linear.

Edge Preserving Properties.

  Median filtering is one kind of smoothing technique, as is linear
  Gaussian filtering. All smoothing techniques are effective at removing
  noise in smooth patches or smooth regions of a signal, but adversely
  affect edges. Often though, at the same time as reducing the noise in
  a signal, it is important to preserve the edges. Edges are of critical
  importance to the visual appearance of images, for example. For small
  to moderate levels of (Gaussian) noise, the median filter is
  demonstrably better than Gaussian blur at removing noise whilst
  preserving edges for a given, fixed window size. However, its
  performance is not that much better than Gaussian blur for high levels
  of noise, whereas, for speckle noise and salt and pepper noise
  (impulsive noise), it is particularly effective. Because of this,
  median filtering is very widely used in digital image processing.

Source: Wikipedia


How Come RIP Guarantees Unique Restoration of the Sparse Solution by $ {L}_{1} $ Minimization?

Compressed sensing theory (see for example Candès & Wakin, 2008) states that, given you have enough measurements $y$, the shown $\ell_1$ minimisation problem recovers $x$ exactly with very high probability. This means that your correct solution must have the smallest possible $\ell_1$ norm.

This is what the figure shows; the "$\ell_1$ ball" intersects the plane given by $\Phi x' = y$ at the point $x'$ with the smallest $\ell_1$ norm. There are, as you mention, other sparse solutions - i.e. the ones where the plane intersects the other two axes. As we can imagine from the figure, these points however have larger $\ell_1$ norm and are therefore not the correct solution.


Kindly suggest what changes I need to make to get a better resolution of eye diagram without increasing the sampling frequency?

There is none, since the resolution of your eye diagram is directly linked to the sampling rate of your observer. In your case, that observation rate is identical to the sampling rate of your signal-producing system.

However, in reality, you'd use an oscilloscope that has a sampling rate much higher than your transmitter.

You can simulate the same – just use a proper interpolator between your modulator's output and your eye diagram's input. Pay attention to the fact that an insufficiently designed interpolation filter will distort your signal – and that is exactly the reality of digital transmitters: the reconstruction filters after the DAC distort the signal, too, and so does the RF channel between a transmitter and the oscilloscope (and of course, the filters and active components inside an oscilloscope do that, too).


pi/4 QPSK modulation/demodulation in GNU Radio

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


why does the additive synthesis method for a triangle wave require amplitude scaling by 8/pi^2?

Are you familiar with how to obtain the Fourier series of an arbitrary periodic function? Because if you are, I'm confused at your confusion.
This is simply the Fourier series expansion of the triangle wave with fundamental frequency $f_0=1/T$. This can be found by taking the inner product of the function with the basis functions. We note that the DC value (the time-average) is zero, and that the function is odd, so it has no cosine terms. So, let's attempt to expand it as the sum of sines or in other words as
$$ x(t) = \sum_{k=0}^{\infty} b_k \sin\left(\frac{2 \pi k t}{T}\right) $$
where
$x(t)$ is the triangle wave defined on the interval $-T/2$ to $T/2$ as
$$ x(t) = \left\{\begin{array}{cc}
\frac{4 t + 2 T}{T} & -\frac{T}{2} < t < -\frac{T}{4} \\
\frac{4 t}{T} & -\frac{T}{4} < t < \frac{T}{4} \\
\frac{4 t - 2 T}{T} & \frac{T}{4} < t < \frac{T}{2}
& \end{array}\right. $$
and elsewhere by the fact that $x(t) = x(t+T)$.
To calculate our Fourier series coefficients, we must compute the integral
$$ b_k = \frac{2}{T} \int_{-\frac{T}{2}}^{\frac{T}{2}} x(t) \sin\left(\frac{2 \pi k t}{T}\right) dt.$$
Because the function we are integrating is a product of two odd functions, it is even, so we can reduce this to
$$ b_k = \frac{4}{T} \int_{0}^{\frac{T}{2}} x(t) \sin\left(\frac{2 \pi k t}{T}\right) dt.$$
We notice that $x(t)$ exhibits even symmetry about the point $\frac{T}{4}$ while even frequency sine waves of the form $\sin(4\pi k t)$ exhibit odd symmetry about the same point. This is quarter wave symmetry and tells us that the even harmonics will all be zero. So, we only need to solve this for odd values of $k$. To that end, we replace $k$ with $2 k + 1$.
We can once again exploit the even symmetry of the odd harmonics to further reduce the range of our integral to
$$ b_{2 k + 1} = \frac{8}{T} \int_{0}^{\frac{T}{4}} \frac{4 t}{T} \sin\left(\frac{2 \pi (2 k + 1) t}{T}\right) dt,$$
and, skipping the calculus (for brevity), we arrive at
$$ b_{2 k + 1} = \frac{8}{\pi^2 (2 k + 1)^2} \cos(\pi k). $$
This reduces to the formula given on Wikipedia. If you don't see it, let me know and I'll be explicit.


Why are anticausal systems even defined?

Of course they don't exist. But we can stop time and use systems that would be non-causal if we hadn't stopped time. Stop time? Yes, just store your data and work offline / non-realtime. Or work on data that are not temporal but something else, for instance spatial. Of course, non-causal is a misnomer in that case, but the term is still used, in analogy with temporal signals and systems.

One common "non-causal" system that is used in offline mode is filtering followed by reversed filtering using IIR filters. This is done in order to eliminate phase distortions introduced by the conventional (causal/realtime) use of the filter. See this answer for the details.


What is the benefit of using symmetric kernel in Fourier transform?

A set of symmetric kernels along Cartesian axes should make the transform consitent with respect to (at least 90 degree) rotations of the input.  If the 2D input yields a particular 2D output, rotating that same 2D input by an angle of 90 degrees, should yield the same 2D output rotated by 90 degrees.  If the kernels in the two dimensions were not symmetric, the transform would yield different output that could not be rotated back to match the original output.

I'd have to think if the invariance applied to angles that were not multiples of 90 degrees.

You may also want to have a look at the Abel transform which has a circularly symmetric kernel, useful for physical situations with circular symmetry.


How to improve the reliability of KF estimation?

For very long times, ive used the KF as a parametrized estimator, training the noise matrices, fitting them standard conditions. Just like a "neural net" on the covariance noises. And you have killer results ;).

And the method?... I've applied any standard nonlinear, gradient descent optimization algorithm, such as Newton Raphson, Levenberg Marquardt...

Sweet, simple, unexpected, and excelent results :). Though unpublishable xd.

If you have data, we can try them...


How to remove ECG artifacts from EMG data?

Here is principal component analysis of the data you supplied. PCA will extract the features within the signal as is done in picture 1.

Note that the signal will be broken down as it is presented in picture 2 from the paper in 3 (page 15). To remove, just note that most of the principal components look like they come from the ecg data, so lets just take all the eigvectors after 10 and project them onto our signal and add them up to get emg data (see picture 4). I'm not entirely sure what the emg data is suppose to look like, but this looks close to diaphragm emg data I see on google images. Try excluding even more eigenvectors to see if you get what you are looking for.

    close all
    clear all
    H = csvread('\example.csv');

%regroup data into snapshots of the signal
%to treat as a group of random processes
l = 2710

Z = (H(1:l*10));
X = reshape(Z,[l, 10]);[![enter image description here][5]][5]
RV = X*X'/length(X)
%eigen decomposition
[U,D,V] = svd(RV);

figure
stem(diag(D))
title('eigenvalues')

%projection
figure
for i  = 1:12
xx0 = V(:,i).*Z(1:l);
subplot(4,3,i)
plot(-xx0)
strr = sprintf('principal component: %d', i);
title(strr)
end

sig = 0;
for i = 11:length(V)
  sig = sig + V(:,i).*Z(1:l);
end

figure
plot(sig)
title('Emg data?')

figure
subplot(121)
plot(V(:,11).*Z(1:l))
title('Principal component 11, emg?')
subplot(122)
plot(V(:,12).*Z(1:l))
title('Principal component 12, emg?')


Why sinc function defined on minus side?

As it relates to linear filtering, the sinc function often comes up because it is impulse response of an ideal lowpass filter (i.e. one that has a perfectly flat passband and a transition region of zero width). In many applications, having such a filter response would be a great tool, as you could perfectly excise any out-of-band interference without disturbing a signal's passband at all.

However, there's a problem, as you noted in your question: $\text{sinc}(t)$ has infinite support! In other words, the impulse response of an ideal lowpass filter is infinitely long. That's not necessarily an issue (a lot of filters have infinitely-long impulse responses), but the big problem is that if the filter is causal (which it certainly must be in order to be realizable), then it will have infinite delay. So, if you were clever enough to build such a contraption, you would find that you would need to wait an infinite amount of time in order to see the filter's response for a given input signal. That tends to make such an idealized filter less useful in practice (even if you could build it).

With that said, while an ideal filter may not be practical, an approximation to an ideal filter can be. One simple way of designing filters, the window method, relies upon taking some filter prototype (an ideal lowpass filter, for instance) and multiplying its impulse response by a finite-length window function. The specific window function used will affect the properties of the resulting filter, but the typical goal is to approximate the prototype filter's response with something that is actually realizable.  As a result, a windowed ideal lowpass filter will still have an impulse response that looks $\text{sinc}$-like over some finite-length interval. Since a causal filter's impulse response starts at $t=0$, the main lobe of the impulse response won't occur until some later time; this is what manifests itself as filter delay.

The plot that you showed is for $\text{sinc}(t)$, which is the ideal, noncausal ideal lowpass filter's impulse response. Translating that to an implementable approximation requires truncating its response to a finite length (i.e. through windowing it) and shifting its response along the time axis so that the filter is causal (after shifting, its response would begin at $t=0$ as you expect).


Deconvolution Question on Article "Deriving Intrinsic Images from Image Sequences" by Yair Weiss

I didn't have time to fully read the article yet I'll try to answer according to your question.

If you take the first equation and convolove it on the left by $ {f}_{i} $ and use the second identity you can see this specific $ g $ holds the equation.

Regarding how to get this $ g $, I think he uses "Deconvolution" in the Fourier domain.


Face recognition using independent component analysis (ICA)

For your data set of images, first vectorize the images by raster scanning them, and making them vectors. Thus, say you have $M$ images, each of size 64*64 pixels. Then the total number of pixels per image is $N=64^2$, which means $N=4096$. Now, you have an image matrix of size $M\text{x}N$.

For this image matrix, what you want to do is find the $M$ independent components, of again, length $N$. In other words, you want to decompose this set of images into another set of images, except those images make up the independent components of the original faces. (Similar to eigen-faces, but different).

To answer your question about the mixture, the set of faces that you start with is the mixture. In other words, the set of faces that you have, is assumed to be a mixture of independent faces that you are trying to find.

ICA will find those independent features for you, and for natural images, independent components turn out to be features with extremely high kurtosis, namely, the edges.


should I apply low-pass filter when calculating central derivative?

No, you are not increasing the sampling time.  The sample rate stays the same.  The Central difference derivative filter is simply longer than the First difference derivative filter.  A longer filter just means more delay before you get your result.

No, you should not low pass filter without a clear reason for why it is needed to meet your signal processing objectives.  Your reasoning is flawed in this case, so don't apply a LPF.

Note that the Ideal derivative filter multiplies a signal's spectrum by $j\omega \cdot e^{j\omega \mu}$ in the frequency domain.  In the time domain, this ideal filter is an infinitely long filter with taps:

$$h[n] = \dfrac{1}{n+\mu}\left[\cos(\pi \left[n+\mu\right])-\mathrm{sinc}(n+\mu)\right]$$

Where

$n$ is the filter tap index in $(-\infty, \infty).$

$\mu$ corresponds to an inter-sample shift of the filter taps in the range $[0.0, 1.0]$ samples.

Given the above expression for $h[n]$, note that:

The Central difference is $h[n]$ truncated to 3 taps with $\mu = 0$.

The First difference is a scaled version of $h[n]$ truncated to 2 taps with $\mu = 0.5$.

Update to add filter magnitude responses plot


you can look at

https://math.stackexchange.com/questions/1290620/what-is-the-importance-of-phase-spectrum-in-fourier-transform

the short answer is that for noise, phase is random.  There shouldn't be any structure to exploit, while for signal, there often is.


By Welch's method, you can calculate the power spectrum by averaging the magnitude of a bunch of FFT frames (or time frames of the STFT). If your signal is $x(t)$, and its STFT is $X(\omega, \tau)$, where $\omega$ is the frequency bin, $\tau$ is the time frame and $T$ is the total number of frames, the PSD is
$$P_x(\omega) = \frac{1}{T} \sum_{\tau=0}^{T-1} |X(\omega, \tau)|^2 $$
When calculating the STFT, you need to make sure you get constant overlap add. This means, with the Hann window, you need a hopSize of N/2 (where N is your window length). I'd also keep the window length shorter, about 256/512 samples. This is not compulsory in your case, because your signal is stationary. Had it not been stationary you should have chosen a smaller window length. Here is your code, modified.
    y, sr = librosa.load('Noise Recorder.wav',sr=None,duration=2)
    n_fft = 8192
    win_len = 512
    S = librosa.stft(y=y,n_fft=n_fft,win_length=win_len,hop_length=win_len/2,window='hann')     # Short FFT
    F = librosa.fft_frequencies(sr=sr, n_fft=n_fft) # Frequency components
    Abs = np.abs(S)  # Same as mag, phase = librosa.magphase(S)

    PSD = np.mean(Abs**2, axis=1).  #this is the power spectrum
    PSD_dB = 10*np.log10(PSD). #power spectrum in dB
    RMS = np.sqrt(np.mean(PSD**2))) #RMS of the power spectrum
    RMS_dB = 20*np.log10(RMS) # RMS of power spectrum in dB

    ##Theoretical RMS (in dB) for a sine wave of amplitude A dB
    A_lin = 10**(A/20.0)           #convert dB to linear amplitude, A = 114 in your case
    RMS_exp = 20*np.log10(A_lin/n_fft)   #expected RMS in dB

`


Advantages and disadvantages in term of radar performance of LFM chirp

Assuming this is not homework, you can summarize LFM performance as

Advantanges

LFM is easy to implement in hardware - even more so if you can use stretch (de-chirp) processing. It has been the most popular type of pulse compression.
Good range resolution which is a function of bandwidth.
Autocorrelation (matched filter output) has good side-lobe levels.
The LFM spectrum is "nice" given its sharp band edges. This make it more robust to signal processing chains with lots of filtering.
LFM is Doppler tolerant. This means that the matched-filter output (think SNR) is robust to Doppler shifts from the return signal.

Disadvantages

LFM exhibits range-Doppler coupling. A Doppler shift will shift the peak of the matched filter output to a different delay (range).


given an image and its fourier transform, find another image fft using fourier properties

Imagine the same problem in 1-dimension and you will be able to find the answer.

Hint: If you know the Fourier transform of a delta function (or time-shifted delta function) and the Fourier transform of a rectangular window, then you should be able to answer the above problem.

P.S. Cannot comment because I need 50+ reputation to do so.


Comparing arithmetic complexity of FFT radix-2 and convolution

Let the input signal $x[n]$ has a legth of $N=50$ samples, and the filter $h[n]$ has a length of $M=10$ samples, then the output $y[n]$ (by linear conv) will have a length of $L = N+M-1=59$ samples.
If you use time-domain convolution, then the number of real MACS can be seen to be  $N \cdot M = 50 \times 10 = 500$.
If you want to use radix-2 FFT to implement the linear convolution result, then you should select a length of $R = 64$ for FFTs. And you will: 1- convert $x[n]$ and $h[n]$ into $X[k]$ and $H[k]$ by two $R$-point FFTS, 2- multiply the results to get $Y[k] = X[k]H[k]$, and 3- apply inverse FFT of $R$-point on $Y[k]$ to get the output $y[n]$.
Each $R$-point FFT (and IFFT) requires about $\frac{1}{2} R  \cdot \log_2(R)$ complex MACs. One complex MAC is 4 real MACs, therefore this is equivalent to $ 2 ~R ~ \log_2(R) = 128 \log_2(64) = 128 \times 6 = 768 $ real MACs. And the total of two FFTs and one inverse FFT requires about $3 \times 768 = 2304$ real MACs. The intermediate multiplication also requires $4 \times R = 256$ real MACs and hence FFT based implementation requires a total of $2560$ real MACs.
So in this case, time-domain convolution is more efficient than an FFT based implementation.


It is usually called DC offset or DC component. Since you use $t$ as a variable, which suggests "time", the quantity $b$ would normally be called frequency. Also note that $c$ is not the phase but a time offset. The phase $\phi$ is usually defined as in

$$x(t)=A\sin(2\pi ft+\phi)\tag{1}$$

Comparing (1) to your formula gives $\phi=-2\pi bc$.


How should I modulate a $16$ bit input using FSK?

Since the comments on the question and the OP's responses
don't seem to be converging, here goes.

Pick a data rate $R$ bits/second (e.g. $10$ kilobits/second)
that is acceptable to your client or boss.
The FSK signal
will then convey one bit every $R^{-1} = T$ seconds (every
$T= 100$ microseconds if you chose the $R =10$ kbps data rate
suggested as an example), and the $16$ bits that you need to
transmit will be transmitted by an FSK signal $x(t)$ that will
have a total duration $16T$ seconds ($1600$ microseconds for the
example).  Specifically, let us assume that you need to define
what $x(t)$ is for all $t$ in the interval $[0, 16T]$.
Number the $16$ bits typed in by the user as $a_0, a_1, \ldots, a_{15}$
where each $a_i$ can have value $0$ or $1$. Include in your
program a snarky error message to be sent to the user if any of
the $a_i$ have value other than $0$ or $1$
Pick two positive integers $m_0$ and $m_1$. The FSK signal
described below will use frequencies $m_0/T$ Hz and $m_1/T$ Hz.
It is convenient, and possibly might even please your client or
boss, if you choose consecutive integers as $m_0$ and $m_1$ such
that the FSK signal fits in the frequency band allocated to the system.
Define the signal $x(t)$ as follows.  For $i = 0, 1, 2, \ldots, 15$,
$$x(t) = A \sin\left(2\pi\left[(1-a_i)m_0 + a_i m_1\left]\frac{t}{T}\right.\right.\right),~~~
iT \leq t < (i+1)T.$$
Several points need to be noticed with respect to this definition:

The quantity in square brackets has value $m_0$ if $a_i = 0$ and value $m_1$ if $a_i = 1$ and so during the time interval $[iT, (i+1)T)$, the signal
consists of either $m_0$ periods or $m_1$ periods of a sinusoid. At the
endpoints of the interval, the signal has value $0$. The signal thus has
the phase continuity that is desirable in an FSK signal as noted
by @JasonR in a comment.
If you are using some programming language to create an array of
$16$ gazillion points that you will then pass to a graphics package to display
as a graph of values of $x(t)$ as a function of $t$, then note that the
formulas used to create the first gazillion points will differ from the
formula used to create the next gazillion points etc.. The
formula needs to use $a_i$ to compute the values of the $i$-th
subset of gazillion points. If you want
to have a nice smooth graph of the sinusoids, it is necessary
to have $m_0$ and $m_1$ be onsiderably smaller than a gazillion
so that you get plenty of points
per period of each sinusoid.  Your graph will also have a long
horizontal axis that will likely extend past one page even if printed in
landscape mode on legal-size paper.


A prototype is an element of the data space that represents a group of elements. On the context of clustering (e.g. under a leaf), a cluster prototype serves to characterize the cluster, their elements. You can have a look at Cluster analysis: basic concepts and algorithms for instance, taken from Introduction to data mining.

In a more specific acceptation, it denotes an instance of clustering for combined numerical and categorical values, in the spirit of k-means, k-modes, k-medians, see for instance Into the world of clustering algorithms: k-means, k-modes and k-prototypes:

  For numerical and categorical data, another extension of these
  algorithms exists, basically combining k-means and k-modes. It is
  called k-prototypes.


This project had fsk demodulation implemented( actually, they reuse another open source project).
I think it would be much easy for you to build python bridge to it, then to recreate all of it from scratch.


Measuring amplitude of a spectral component

A variation on @geometrikal's answer would be to use the Goertzel Algorithm.

$$
y[n] = x[n] + e^{+j\omega_0} y[n-1]
$$
where $x[n]$ is your input signal, $\omega_0$ is the frequency you know, and $y[n]$ is the (complex-valued) estimate of the discrete Fourier coefficient.


OK, I solve the problem. The problem was that I was overriding the files when I used the "sox" command, as I want them to convert into RIFF format. However, "sox" didn't convert them correctly when using the same file name, i.e. when overriding the files. When I changed the file name of the overridden files, it converted them into RIFF format.


Index of stationarity of a time domain signal

This a very complicated question, and I would say a still open topic. The concept of stationarity is manifold, from pure statistics to applied DSP (strict, strong, wide-sense, quasi-stationarity, cyclo-stationarity, to refer to a recently closed question). The lack of access to sound models and faithful realizations renders the quest quite difficult.

Non-stationary is a non-property (like non-linearity, non-causality). Refocusing on some breeds may help its characterization, whether  you are interested in:

randomly stepping piece-wise constant signals
"sines" of  varying frequency (like chirps)
stochastic data of changing distribution
etc.

Yet, a common framework resides in using time-frequency representations, assuming some local stationarity and trying to assess the non-stationarity according to some metric/distance/divergence (norms, Kolmogorov, Kullback, Jensen, etc.) over the time. Here are two possible (first) references:

Stationary index for abrupt changes detection in the time-frequency plane, 1998
Testing Stationarity With Surrogates: A Time-Frequency Approach, 2010


Part of the advantages of using dB to express gain or attenuation is that they're inherently logarithmic.
So if you're speaking in dB, when you add a 10dB attenuation to another 10dB of attenuation, what you're saying is that you're multiplying the signal by $\sqrt{0.1}$, and then by $\sqrt{0.1}$ again.  The result is 20dB of attenuation -- or multiplying the signal by $\sqrt{0.01} = 0.1$.
So if your ultimate gain is $A$, as a normal number, and you have attenuations $a_1, a_2, \cdots$ in dB, then the actual $A$ by which you multiply your signal is $$A = 10^{-\frac{1}{20}\sum a_n}$$


The actual "data" in the IQ stream from the USRP would depend on what frequency and signal you were looking at when the IQ stream was captured and presumably saved to your laptop.

This picture shows that you can filter the I and Q completely separately - so rather than thinking of them as complex values think of them of as two streams of real numbers.

If your goal is specifically to filter them using a hamming window (which is typically a frequency domain concept) then taking the FFT of the data is appropriate.

When first introduced to frequency domain transforms the example are almost always real-valued signals. However, it is very appropriate (and even more useful) to take FFTs of complex valued signals. In this case - leaving them as "complex" symbols and using an FFT library to take the FFT would be helpful. MATLAB, SciPy for Python, and Octave all have FFT implementations and several of them are based around the FFTW library.

To "filter" your IQ samples with the hamming window you'll probably want to multiply the FFT of the IQ point-by-point with the Hamming Window. In MATLAB this would look something like fft(IQData).*HammingWindow - where hamming window is a purely real vector and the FFT(iq) is a complex stream.

Your final question regarding finding the frequency content is probably best visualized by taking the magnitude/norm/abs of the IQ samples. Taking sqrt(I^2+Q^2) will take the two complex numbers and turn them into a single real value that can be plotted vs. frequency. If you want to know how to scale your x-axis of this plot (assuming y-axis is the sqrt(I^2+Q^2)) then this link should be helpful (and short answer is matlab's fftshift() function helps).

http://www.bitweenie.com/listings/fft-frequency-axis/

-Arch-


OFDM signal spectrum and characterisitics

Your statement isn't strictly correct; the spectrum of an OFDM signal isn't a rectangular pulse. Think about it this way: OFDM is really just a method to efficiently modulate multiple low-rate narrowband signals that are sent simultaneously. Therefore, the OFDM signal's spectrum is really just the sum of the spectra of each of its narrowband components.

For example, if your OFDM signal has 64 BPSK subcarriers, then its spectrum is equal to the spectra of 64 narrowband BPSK signals, spaced at the OFDM subcarrier spacing. If you were to estimate the OFDM signal's power spectrum with a wide enough resolution bandwidth, then the spectrum would appear to be rectangular. However, if you look closely enough at the signal's spectrum, you'll notice that it isn't flat.

Here is an example MATLAB script that demonstrates the phenomenon:

% generate random BPSK symbol vector; this is equal to the "low-resolution"
% spectrum, calculated with bin spacing synchronized with the OFDM
% subcarrier grid
X = 1 - 2 * (randn(64,1) > 0.5);
f = -0.5 + (0:length(X)-1) / 64;
% create time-domain OFDM symbol using the IFFT
x = ifft(X);
% calculate interpolated spectrum to see what it looks like between OFDM
% subcarrier centers
Sh = fftshift(fft(x,4096));
fh = -0.5 + (0:length(Sh)-1) / 4096;

figure;
plot(f,abs(X)); hold all; grid on;
plot(fh,abs(Sh));
legend('Regular', 'Interpolated');

As the modulated OFDM signal doesn't have a perfectly rectangular spectrum as you guessed, it does not necessarily have a sinc-like shape in the time domain either.


Zero-crossing rate (ZCR) might be useful for voiced/unvoiced frame discrimination, speech/music discrimination, but it is of much lesser importance in speech recognition. One reason is that it is pitch-dependent and not robust to background noise or hum. It is not difficult to craft very different signals (say a female voice saying a phone and a male voice saying another phone) that have the exact same ZCR; or to add some noise to a signal in a way that will half or double its ZCR.

Using zero-crossing rate might work in the single-speaker + prior training case. Collect a handful of samples of each words, compute the zero crossing rate, and train a naive bayesian classifier for this two-class problem. Be prepared to poor performances, and don't expect your system to generalize to other speakers.

The only plausible reason one would use such a feature would be on embedded systems with very limited computational power. Otherwise, there's no reason not to use "proven" features like MFCC or AR coefficients.


Best digital noise reduction in optical spectra measurement

For the additive white noise channel, it can be shown that matched filtering maximizes SNR; and for constant illumination, the "transmitting pulse filter" is a rectangle, and the match to that is a rectangular weighting, i.e. an average, too.

The idea is simply that you "my signal is a vector; multiplication / convolution then becomes a dot product. Cauchy-Schwarz Inequality gives me an upper bound for the norm of signal·receive window / filter, and I reach that when what I do at the receiving end is a linear multiple of what I do at the transmitting end".

For unknown spectral distributions, you'd try to find a filter on the receiving end of things that is as orthogonal as possible to the spectrum of the noise. You'd need to estimate that first, though! Devise some calibration method; but, really, when you're observing the same spot for long enough, maybe just look at the spectrum and ignore the constant value?

Assuming you have some non-white noise, you can "abuse" techniques from the world of spread-spectrum communications design. Imagine the following:

You can modulate the intensity of your light source to two values, a and b, over time. Imagine sending a pseudo-random sequence consisting of times when you send the higher intensity, and of the lower intensity:

a      ----    --------    ----

(a+b)/2

b          ----        ----

At the receiver, the first thing you do is simply subtract the average of the received signal intensity; that way, what you receive becomes:

c      ----    --------    ----

0

-c         ----        ----

Now, that's a very simple spread spectrum signal (in fact, DSSS), and it's almost certain that its spectral properties are pretty different from these of noise.

Now, you just correlate with the high-low sequence that you've sent, considering "high" a "1", and "low" a "-1". When you shift that "+1 -1 +1 +1 -1 +1" sequence correctly, you'll see that by point-wise multiplication between "c -c +c +c -c +c" (which you received), you get a very high amplitude: 1·c + (-1)·(-c) + 1·c + 1·c + (-1)·(c) + 1·c = 6c . In reality, you use lengths much higher than 6!

You still get the same gain as when averaging white noise over a length of $N=6$, but for low-frequency-heavy noise distributions, you'd get noise suppression, too. (you buy that by getting less noise suppression at higher frequencies! No such thing as a free lunch; but you might want to do that, because you might have other light sources in your room – for example, 100 Hz modulating room light, which you can eliminate efficiently).

A bit of general commentary:

That averaging gives you the gain you state depends on the assumption that noise is uncorrelated with signal! In general, that's not the case for you: The energy of noise is very much dependent on current forward current in a semiconductor junction; that's kind of a bummer. What the spread-spectrum approach does is in fact a bit of noise shaping, and it might well suppress a bit of the harmonics that nonlinearity introduces, but I'm not expert enough in photonics to tell you how severe these nonlinearities are.

Also, I don't know the reflectivity of what you're measuring. Just be warned that white LEDs aren't really white: they consist of multiple light sources that generate different wavelengths; same goes for your photo sensors: They aren't uniformly sensitive to all wavelengths. If you haven't already, calibrate your sensor with a color proof sheet!


Common Use Cases for 2D Non Separable Convolution Filters

There are many.  In general, steerable filters are not separable.  Few are, but many are not.

Here is a short list:

Gabor filters in certain directions.  Of course, one can use the bases that are separable due to the fact that Gabor filters are steerable.  Gabor filters have the advantage that they are optimally localized in time and frequency.  As a degenerate case, consider a non-symmetric 2D Gaussian where the principal components are not horizontal or vertical.
Difference of gaussian (DoG) used to model retinal ganglion and thalamic receptive fields are also not separable.
Marr-Hildreth for edge detection is not separable.
A bunch of filters used for segmentation are not separable: http://dl.acm.org/citation.cfm?id=720464&CFID=260338827&CFTOKEN=99199893

The list can go on and on, and these are just common examples!  You can certainly design FIR filters to fit your need, and there are methods to do this.  That is the whole purpose of a large sub-field of DSP.

In general when you design a 2-D filter, if it is anisotropic and not oriented horizontally or vertically, chances are that it is not separable.  For example, if you were to design a low pass (pre-)filter to remove high frequency that may cause aliasing in a hexagonal sampling system.  This would have to be an analog filter, but for digital, consider a down sampling scenario.  For optimal design, you'd want to keep as much information as possible, so your filter would most likely match the hexagonal sampling grid.  This will definitely give you a filter that is not separable.


Why do we say that "zero-padding doesn't really increase frequency resolution"

Resolution has a very specific definition in this context. It refers to your ability to resolve two separate tones at nearby frequencies. You have increased the sample rate of your spectrum estimate, but you haven't gained any ability to discriminate between two tones that might be at, for instance, 236 Hz and 237 Hz. Instead, they will "melt together" into a single blob, no matter how much zero-padding you apply.

The solution to increasing resolution is to observe the signal for a longer time period, then use a larger DFT. This will result in main lobes whose width are inversely proportional to the DFT size, so if you observe for long enough, you can actually resolve the frequencies of multiple tones that are nearby one another.

--

To see how this plays out, here's a plot of the zoomed-in FFT of the addition of two signals: your original sinusoid, and one that differs in frequency from it by 0 to 100 Hz.

It's only towards the 100Hz difference end of the plot (left-hand side here) that you can distinguish (resolve) the two.

Scilab code for generating the plot below.

f = 236.4;
d = 10;
N=441;
fs=44100;
extra_padding = 10000;

t=[0:1/fs:(d/1000-1/fs)]
ff = [0:(N+extra_padding-1)]*fs/(N+extra_padding);

x = sin(2*%pi*f*t);

XX = [];

for delta_f = [0:100];
    y = sin(2*%pi*(f+delta_f)*t);
    FFTX = abs(fft([x+y zeros(1,extra_padding)]));
    XX = [XX; FFTX];
end

mtlb_axis([0 1300 0 500])

figure(1);
clf
[XXX,YYY] = meshgrid(ff,0:100);
mesh(XXX(1:100,[50:90]),YYY(1:100,[50:90]),XX(1:100,[50:90]))


How to simulate High Power Amplifier(HPA), like Solid State Power Amplifier (SSPA)

The model works fine.  You need to slow down and check your (or someone else's? ) steps.  I just worked through one pass of this, without pskmod, and only one value was high enough voltage to push the model into a nonlinear region.  Work through it step by step, you'll see that your x values are too low in magnitude to really get out of the linear region.  The x values have to approach A to get distortion.


Filter Coefficients in MATLAB are always float data type

Normally you would use fixed-point arithmetic for performing DSP on a a floating-point-challenged CPU.

So for, say, 16 bit fixed-point coefficients with a range of -1.0 to +1.0 your coeffs would translate to e.g.:

0.465    =>    0.465 * 0x7fff = 0x3b5

If your microcontroller does not have direct support for fixed point arithemtic (most DSPs do, general purpose microcontrollers typically do not) then you'll need to take care of scaling when you multiply or divide. (Addition and subtraction work as normal of course, provided you are not mixing different fixed point types.)


White noise has no component at any specific frequency since it is a noise density and therefore has energy across a range of frequencies. The range of any given specific frequency, including DC is zero, and therefore the energy at DC (and any specific frequency) is indeed zero.

In practice you cannot actually observe DC or any specific frequency as a single point as that would require an infinite observation time. However, the longer you observe, the narrower will be the effective range of frequencies you observe and for this reason the mean will be non-zero but approach zero as your observation time approaches infinity.


Is the Canny algorithm enough for creating a feature descriptor image for input to an SVM?

Yes and no.

Yes, the binary image output from the Canny edge detector is a feature set.

As the Feature (computer vision) wiki page describes:

  ... an edge can be represented as a boolean variable in each image point that describes whether an edge is present at that point. Alternatively, we can instead use a representation which provides a certainty measure instead of a boolean statement of the edge's existence and combine this with information about the orientation of the edge.

What you have is the former; a boolean variable for every pixel in the image, true where there is an edge. A Hough Transform would be an example of a method to get the latter; edges represented by measure and orientation.

If you take the same image and convert it to edges, the feature set will be the same. If you take two different images and convert them to edges, the feature set will likely be different. So yes, you can model an image using the Canny edge detector alone.

Now, no, the binary image output from the Canny edge detector is probably not a very useful feature set.

Consider a single training image of an apple and three test images: an apple, an orange, and a banana. By shape, what we are focusing on by focusing on the Canny edge detector, we would expect the test apple to be most closely related to the training apple, the orange second, and the banana last.

In a candy land of image analysis, where the images are all perfectly composed, their edges could line up such that the apples edges overlap exactly, the orange's edge overlaps a lot with the apple, and the banana's edges only overlap for the few pixels they intersect. However, in reality you are probably comparing images of different apples with different backgrounds and different composition. Imagine the test apple is shifted one pixel down and to the right, now it hardly overlaps with the training apple at all; the orange will be a better match. In this example, the binary image output from the Canny edge detector is not a robust feature descriptor.

Think of the Canny edge detector algorithm as a transformation of the image. It can be a useful step in getting to a powerful set of features which describe your image, but it itself is not usually a powerful feature descriptor.

So, what would be a powerful feature descriptor?

..well you didn't ask that ;)


I read that the coefficients are always 1 for moving average so this is what I did - is that correct?

Technically no. For a moving average filter of length $N$ all the coefficients are $1/N$. In practice most people use $1$ for all coefficients and then scale the output with $1/N$ It gives the same result and it's typically more efficient, but if you were to look at the difference equation, the coefficients would show up as $1/N$

Moving average filters are symmetrical,

Depends a bit exactly how exactly you define symmetry. You can certainly time reverse it without changing anything.

have no L-th band filter,

I have no idea what you are mean by "L-th band filter"

and are single rate filters right?

The concepts of multi-rate and moving average are mostly orthogonal. You can use the together , individually or not all.


Static gain refers to the DC gain. Namely, it would be the ratio of the output and the input under steady state condition.

Due to DC corresponding to $\omega=0$, in the $z$-domain DC would correspond to $z=1$ because $z=re^{j\omega}$, with $r=1$ (i.e. you have to evaluate your transfer function $G(z)$ in the unit circle to get the frequency response).

EDIT:

The frequency response of a system is represented in the Fourier domain, $G(\omega)$. The DTFT of a given discrete sequence can be calculated as:

$$G_F(\omega)=\sum_{n=-\infty}^{\infty}g[n]e^{-j\omega n}$$

On the other hand, the $z$-transform is defined as:

$$G_Z(z)=\sum_{n=-\infty}^{\infty}g[n]z^{-n}$$

It's easy to see that if we evaluate $z=e^{j\omega}$ (the unit circle), then the $z$-transform returns the frequency response.


Additivity requires a little more than a direct addition of some $x_1$ to some $x_2$. It should involve any linear combination of  $x_1$ and $x_2$, $a_1 x_1 + a_2 x_2$ where $a_1$ and $a_2$ are "kind of scalars" (this could make us dive into complicated stuff, like in  Does scaling property imply superposition?)

In your case, however, linearity works, because you are touching the "time base" ($t$ vs. $e^t$) without touching the amplitudes. Your operation is a sort of time warping.

Accordingly:

$$\mathcal{S}(a_1 x_1(t) + a_2 x_2(t)) = a_1 x_1(e^t) + a_2 x_2(e^t) = \mathcal{S}(a_1 x_1(t)) + \mathcal{S}(a_2 x_2(t))  $$


What Is the Relationship Between a Kalman Filter and Polynomial Regression?

1. There is a Difference in terms of optimality criteria
Kalman filter is a Linear estimator. It is a linear optimal estimator - i.e. infers model parameters of interest from indirect, inaccurate and uncertain observations.
But optimal in what sense? If all noise is Gaussian, the Kalman filter minimizes the mean square error of the estimated parameters. This means, that when underlying noise is NOT Gaussian the promise no longer holds. In case of nonlinear dynamics, it is well-known that the problem of state estimation becomes difficult. In this context, no filtering scheme clearly outperforms all other strategies. In such case, Non-linear estimators may be better if they can better model the system with additional information. [See Ref 1-2]
Polynomial regression is a form of linear regression in which the relationship between the independent variable x and the dependent variable y is modeled as an nth order polynomial.
$$ Y = a_0 + a_1x + a_2x^2 + \epsilon $$
Note that, while polynomial regression fits a nonlinear model to the data, these models are all linear from the point of view of estimation, since the regression function is linear in terms of the unknown parameters $a_0, a_1, a_2$. If we treat $x, x^2$ as different variables, polynomial regression can also be treated as multiple linear regression.
Polynomial regression models are usually fit using the method of least squares. In the least squares method also, we minimize the mean squared error. The least-squares method minimizes the variance of the unbiased estimators of the coefficients, under the conditions of the Gauss–Markov theorem. This theorem, states that ordinary least squares (OLS) or linear least squares is the Best Linear Unbaised Estimator (BLUE) under following conditions:
a. when errors have expectation  zero i.e. $E(e_i) = 0 $
b. have equal variances i.e.  $ Variance(e_i) = \sigma^2 < \infty $
c. and errors are uncorrelated i.e. $ cov(e_i,e_j) = 0 $

NOTE: that here, errors don't have to be Gaussian nor need to be
IID. It only needs to be uncorrelated.

2. Kalman Filter is an evolution of estimators from least square
In 1970, H. W. Sorenson published an IEEE Spectrum article titled "Least-squares estimation: from Gauss to Kalman." [See Ref 3.] This is a seminal paper that provides great insight about how Gauss' original idea of least squares to today's modern estimators like Kalman.
Gauss' work not only introduced the least square framework but it was actually one of the earliest work that used a probabilistic view. While least squares evolved in the form of various regression methods, there was another critical work that brought filter theory to be used as an estimator.
The theory of ﬁltering to be used for stationary time series estimation was constructed by Norbert Wiener during 1940s (during WW-II) and published in 1949 which is now known as Wiener filter. The work was done much earlier, but was classiﬁed until well after World War II). The discrete-time equivalent of Wiener's work was derived independently by Kolmogorov and published in 1941. Hence the theory is often called the Wiener-Kolmogorov filtering theory.
Traditionally filters are designed for the desired frequency response. However, in case of Wiener filter, it reduces the amount of noise present in a signal by comparison with an estimation of the desired noiseless signal. Weiner filter is actually an estimator. In an important paper, however, Levinson (1947) [See Ref 6] showed that in discrete time, the entire theory could be reduced to least squares and so was mathematically very simple. See Ref 4
Thus, we can see that Weiner's work gave a new approach for estimation problem;  an evolution from using least squares to another well-established filter theory.
However, the critical limitation is that Wiener filter assumes the inputs are stationary.  We can say that Kalman filter is a next step in the evolution which drops the stationary criteria. In Kalman filter, state space model can dynamically be adapted to deal with non-stationary nature of signal or system.
The Kalman filters are based on linear dynamic systems in discrete time domain. Hence it is capable of dealing with potentially time varying signal as opposed to Wiener. As the Sorenson's paper draws parallel between Gauss' least squares and Kalman filter as

...therefore, one sees that the basic assumption of Gauss and Kalman
are identical except that later allows the state to change from one
time to next. The difference introduces a non-trivial modification to
Gauss' problem but one that can be treated within the least squares
framework.

3. They are same as far as causality direction of prediction is concerned; besides implementation efficiency
Sometimes it is perceived that Kalman filter is used for prediction of future events based on past data where as regression or least squares does smoothing within end to end points. This is not really true. Readers should note that both the estimators (and almost all estimators you can think of) can do either job. You can apply Kalman filter to apply Kalman smoothing.
Similarly, regression based models can also be used for prediction. Given the training vector, $X_t$ and you applied $Y_t$ and discovered the model parameters $α_0 ... a_K$ now for another sample $X_k$ we can extrapolate $Y_K$ based on the model.
Hence, both methods can be used in the form of smoothing or fitting (non-causal) as well as for future predictions (causal case). However, the critical difference is the implementation which is significant. In case of polynomial regression - with entire process needs to get repeated and hence, while it may be possible to implement causal estimation but it might be computationally expensive. [While, I am sure there must be some research by now to make things iterative].
On the other hand, Kalman filter is inherently recursive. Hence, using it for prediction for future only using on past data will be very efficient.
Here is another good presentation that compares several methods: Ref 5
References

Best Introduction to Kalman Filter - Dan Simon Kalman Filtering Embedded Systems Programming JUNE 2001 page 72

Presentation: Lindsay Kleeman Understanding and Applying Kalman Filtering

H. W. Sorenson Least-squares estimation: from Gauss to Kalman IEEE Spectrum, July 1970. pp 63-68.

Lecture Note MIT Course ware - Inference from Data and Models (12.864) - Wiener and Kalman Filters

Presentation Simo Särkkä From Linear Regression to Kalman
Filter and Beyond Helsinki University of Technology

Levinson, N. (1947). "The Wiener RMS error criterion in filter design and prediction." J. Math. Phys., v. 25, pp. 261–278.


MATLAB - How do I get the Filter Function to get the output Energy Spectrum?

What u plotted is the complex fft samples.  In order to see the psd in matlab use plot([-1/2:1/ffteln:1/2-1/fftlen]*fs, 10*log10( abs( fftout)./length(fftout)).^2)).  I am assuming of coarse that your signal is real and therefore plotting from -fs/2:fs/2 is standard practice.  For complex signals usually you want to plot from 0:Fs because nyquist is Fs for complex.


Advancing a causal IIR filter's impulse response

Let's consider a causal IIR filter with the following transfer function:

$$H(z)=\frac{B(z)}{A(z)}=\frac{\displaystyle\sum_{n=0}^{M}b[n]z^{-n}}{1+\displaystyle\sum_{n=1}^{N}a[n]z^{-n}}\tag{1}$$

Note that $M$ need not be equal to $N$. Let $K$ be the number of samples that must be removed from the impulse response corresponding to the system $(1)$, so the new desired impulse response is

$$\tilde{h}[n]=h[n+K],\qquad n=0,1,\ldots\tag{2}$$

Let $H_K(z)$ denote the FIR transfer function corresponding to the first $K$ samples of $h[n]$:

$$H_K(z)=h[0]+h[1]z^{-1}+\ldots +h[K-1]z^{-(K-1)}\tag{3}$$

The transfer function of the new filter is then given by

$$\tilde{H}(z)=z^K\left[\frac{B(z)}{A(z)}-H_K(z)\right]=\frac{z^K\big[B(z)-H_K(z)A(z)\big]}{A(z)}\tag{4}$$

As expected, the denominator polynomial remains unchanged, we just get a new numerator polynomial

$$\tilde{B}(z)=z^K\big[B(z)-H_K(z)A(z)\big]\tag{5}$$

Let's first consider the case $K=1$. In this case, the FIR transfer function $H_1(z)$ defined in $(3)$ is trivially given by

$$H_1(z)=h[0]=b[0]\tag{6}$$

Consequently, the new numerator polynomial becomes

$$\tilde{B}(z)=\big(b[1]-b[0]a[1]\big)+\big(b[2]-b[0]a[2]\big)z^{-1}+\ldots +\\+\big(b[N]-b[0]a[N]\big)z^{-(N-1)}\tag{7}$$

if $M=N$. If $M>N$, there are additional terms $b[k]z^{-(k-1)}$, $k=N+1,\ldots,M$, and if $M<N$, there are additional terms $-b[0]a[k]z^{-(k-1)}$, $k=M+1,\ldots,N$.

The computation of the new numerator coefficients according to $(7)$ is very simple and efficient. For $K>1$ this procedure can be applied iteratively. In Matlab/Octave this is (almost) a one-liner:

for k = 1:K
    b = [b(2:M+1);zeros(N-M,1)] - b(1)*[a(2:N+1);zeros(M-N,1)];
    M = length(b) - 1;
end

Note that for any $K$, the order of the new numerator polynomial $\tilde{B}(z)$ equals $\max(M-K,N-1)$.

Original Answer (correct, but the presentation and solution given above is probably better; the examples below are still valid and hopefully useful)

A solution can be obtained by recomputing the numerator coefficients as follows:

$$\tilde{b}[n]=(\tilde{h}\star a)[n],\qquad n=0,1,\ldots,L-1\tag{3}$$

where $\star$ denotes convolution, $L=\max(M,N)$, and $\tilde{h}[n]$ is related to the impulse response $h[n]$ of the original IIR filter by

$$\tilde{h}[n]=h[n+1],\qquad n=0,1,\ldots,L-1\tag{4}$$

The method can be summarized as follows:

Compute $L=\max(M,N)$ coefficients of the original impulse response $h[n]$ for indices $n=1,2,\ldots,L$
Compute the new numerator coefficients from the impulse response and from the denominator coefficients using Eq. $(4)$

I chose a random filter to show the result of the suggested procedure in the figure below (blue: original impulse repsonse: red: new impulse response shifted to the right by one sample for comparison):

The left column below shows the original numerator coefficients alongside the modified numerator coefficients, aligned in such a way that it becomes obvious that for higher indices the new coefficients equal the original coefficients. (This is the case because in this example $M>N$; c.f. Eq. $(7)$ of the new answer above, including the text below Eq. $(7)$).

  -0.67010
   0.65673   0.92477
  -0.61219  -0.74621
  -0.60645  -0.53944
  -0.50316  -0.50316
   0.82855   0.82855
   0.20396   0.20396
   1.15680   1.15680
  -0.29830  -0.29830
   0.10605   0.10605

The denominator coefficients of both filters are

a = [1, 0.4, -0.2, 0.1]

EDIT:

This method can be generalized to remove $K$ initial values of the original impulse response, even if $K>M$:

Compute $L=\max(M,N)$ coefficients of the original impulse response $h[n]$ for indices $n=K,K+1,\ldots,K+L-1$
Compute the new numerator coefficients according to Eq. $(4)$ with $\tilde{h}[n]=h[n+K],\qquad n=0,1,\ldots,L-1$

The resulting numerator coefficients will generally have several trailing zeros (or values sufficiently close to zero) which can be removed.

The following example shows the result of removing $20$ initial values of the impulse response of the filter with coefficients

[b,a] =

   1.000000   1.000000
   2.000000  -0.187178
   3.000000  -0.075960
  -1.000000   0.812250

The numerator coefficients of the new filter are

   0.95236
   0.90498
   0.42252

Of course, these coefficients can no longer be interpreted as (an approximation of) a shifted version of the original numerator coefficients.


An ideal sinusoid does indeed have zero bandwidth, but it is not useful to transmit one because it also imparts zero information to the receiver.  Once the receiver knows what the sinusoid phase is it can perfectly predict the rest of the signal- thus, no more information is received.

For information to be transmitted there has to be an element of "not knowing".  For digital communications that generally means not knowing if the next bit will be a 1 or 0.  I will try to illustrate where the bandwidth comes from using BPSK, but the exact same principles apply to QAM.

The BPSK example above is ideal in the sense that it is very easy for the receiver to recover.  The receiver can sample the waveform anywhere within a symbol ("symbol" being the waveform representation of a digital bit) and get the full symbol power.  You can model the example BPSK waveform above as a bunch of individual square waves of varying widths and locations (in time) added together.  If you recall your Fourier theory, the Fourier transform of a square is the sinc function.  Thus, in the frequency domain, the total signal is a bunch of sinc functions of varying widths added together.  That makes the resulting waveform's frequency spectrum look a lot like a sinc function (note that this is not 100% true since there is some destructive interference that occurs, but for the sake of this discussion let's ignore that).

Sinc functions, unfortunately, aren't too great in terms of bandwidth.  The figure below shows the absolute value of a sinc function.

This doesn't look too bad until you look at the log scale of the energy, which is shown below.

As you can see from the plot above, at the location of the signal's potential nearest neighbor, the sinc function is only about 13 dB down.  The rolloff from there is quite slow, with it only being about 21 dB down at the fifth side lobe.  This introduces a lot of interference into neighboring signals, and when your receiver filters out the neighboring signals it is also filtering out a lot of your transmitted power, so that power is wasted.

Transmitters employ pulse shaping that rounds out the harsh transitions of the square waves and makes the BPSK waveform look more like a sinusoid.  The following is the same BPSK waveform but using raised-cosine pulse shaping.

And this is the energy vs. frequency view of the raised-cosine pulse.

The main lobe of the raised-cosine and the square wave are essentially the same.  The difference is in the rolloff of the side lobes.

So long story short (too late!) the bandwidth comes from the "random" transmission of ones and zeros.  The minimum bandwidth is determined by the symbol period (shorter period is faster, thus wider bandwidth), and the sidelobes outside of that main bandwidth are minimized via pulse shaping.


Phase response of a FIR Hilbert transformer

The slope represents a delay of $N$ samples. You can center the impulse response to time zero by multiplying the transfer function by $e^{i\omega N}$ (transfer function of a negative delay of $N$ samples) or np.exp(1j*w*N) in Python:

from scipy import signal
import numpy as np
N = 10
b = [2*np.sin(np.pi*n/2)**2/(np.pi*n) for n in range(1, N+1)]
b = np.append(-1*np.flip(b, 0), np.append([0], b))*np.hamming(2*N+1)
w, h = signal.freqz(b)
h = h*np.exp(1j*w*N)
import matplotlib.pyplot as plt
fig = plt.figure()
plt.title('Digital filter frequency response')
ax1 = fig.add_subplot(111)
plt.plot(w, 20 * np.log10(abs(h)), 'b')
plt.ylabel('Amplitude [dB]', color='b')
plt.xlabel('Frequency [rad/sample]')
ax2 = ax1.twinx()
angles = np.angle(h)
plt.plot(w, angles, 'g')
plt.ylabel('Angle (radians)', color='g')
plt.grid()
plt.axis('tight')
plt.show()

Most of the above is from Scipy documentation which you link to. This gets plotted:

Figure 1. Frequency response of the Hilbert transformer centered at time zero.

The phase goes to zero at frequency zero, because the phase frequency response transitions from $-\pi/2$ for positive frequencies to $\pi/2$ for negative frequencies.


Note that $1536 = 3\cdot 512$. You can use a decimation in time method. See http://www.altera.com/literature/an/an480.pdf

I'll add some extra information here in response to your comment to sansuiso's answer (even though this comment actually referred to the pdf file I mentioned above in this answer). As explained in the above document, you need to compute 3 FFTs of length 512. You do this by taking time samples 0,3,6,9, ... and compute their FFT, then you take the samples 1,4,7,10,... and compute another FFT. Finally you compute the FFT of samples 2,5,8,11,... Then you combine the 3 FFTs using the appropriate twiddle factors to obtain the desired 1536-point FFT. This last step is the radix-3 stage. Have a look at Equation (2) on page 3 of the document (also Fig. 1 on page 4).

Also check out this blog post. It also includes a little Matlab program.


Using neural networks for rotation invariant OCR

Aha, this is a quite funny story that I've heard about CV. Are you a bio guy? Any way, here are my suggestions.

If you are a bio guy and just want to finish this project ( I mean successfully identify each insect in a video frame ), go for barcode, QR etc. They are labels though their contents are not directly readable by eyes. However, you will have a much better chance to handle all kind of problems you might encountered when you use manual labels, e.g. occlusion, error correction, lighting etc. There are libraries can decode these labels, and you simply call them.

In this way, you will save a lot of energy to not working on something you are in lack of experience, but use the state-of-the-art technologies.

If you are not a bio guy or these manual labels are not something can be changed, you should first locate these labels and then recognize them. Donot mess up these two things. The first task is detection, and the second one is recognition. In other words, you need to do two things:
a). Train a text/none-text classifier. In this classifier you can use those rotational invariant features.
b). Train a character classifier ( e.g. using NN ) to decide the class of a character candidate you found from (a).

By the way, I donot believe you need to train a NN accepting character at an arbitrary angle. Instead, all you need is to train a NN accepting normal characters, and normalize a candidate using image moments or other techniques before you feed it into the classifier. Or simply rotate your candidate to a certain degree (e.g. for every 10 degree) and then feed them into your classifier and pick the result with the highest confidence.


Depends mostly on your playback headset. If it's diffuse field equalized, use the diffuse field set. For most headsets, this will be the best choice.

The only difference is overall equalization, so you can turn one into the other with a simple filter that gets applied to all HRTFs


How to Mesure the smoothness of a signal

Hurst exponent is a good option.

For a better answer you must better explain your objective. For example if you are comparing signals and if they have the save length, amplitude, sampling frequency and/or others.

If you are looking for very simple approaches then, for example, you can:

calculate the discrete derivative
calculate absolute value
sum all values

Another option could be:

apply a low pass filter
calculate difference between original signal and filtered signal
a smaller difference means a smoother signal.


If

$$h(t)=\frac{1}{1+t^2}\tag{1}$$

denotes the filter's impulse response, then the filter is an infinite impulse response filter, simply because $h(t)$ has infinite support, i.e. it doesn't vanish anywhere (other than for $t\rightarrow \infty$). Note that if the definition $(1)$ of $h(t)$ is valid for all $t$, then the filter is not causal because $h(t)=0$ for $t<0$ doesn't hold.

Since the impulse response is real-valued and it satisfies $h(t)=h(-t)$, its Fourier transform must be real-valued (and symmetrical). Since the Fourier transform is real-valued, the filter's phase is zero.

If you sample that impulse response using $5$ taps at times $t_{-2}=-2T$, $t_{-1}=-T$, $t_0=0$, $t_1=T$, $t_2=2T$, you get an FIR filter, but the zero-phase property remains unchanged due to the symmetry of the filter coefficients. For implementation, the filter impulse response needs to be shifted by two samples to make it causal. This causes a delay of two samples, which changes the filter's phase $\phi(\omega)$ from zero to a linear phase $\phi(\omega)=-2\omega$.


Relation between frequency and phase in PLL?

no, frequency and phase are not independent quantities.  phase is the integral of frequency.  and in the VCO (or NCO) inside a PLL, there is an integrator.  the frequency is proportional to the input of the VCO or NCO, but the phase (which is what goes into the phase discriminator) is the integral of that.

with feedback, a PLL is a servo-control mechanism.  but because of the inherent integrator in the loop, the "controller" is either a PI or an I.  not much of a D in that controller inside a PLL.


Equalization knowing the Channel State Information

To 1.
From the channel impulse response you can build the channel matrix and then compute the ZF equalizer coefficients using the pseudo inverse of the channel matrix.

Once you have the ZF coefficients you can convolve them with the received signal.

To 2.
I guess you are referring here to a PNC system. The best solution is to use OFDM to solve this problem. In this case, you can use Maximum Likelihood detection per subcarrier by checking which combination minimises the Euclidean distance.

To 3.
You can find the equations for ZF equalizer in the book by J. Kurzweil, An Introduction to Digital Communications, 1st Edition.

BTW, the performance of the time-domain ZF equalizer is not very good in frequency selective channels. You are better off with the decision feedback
equalizer.


What is the scaling function and wavelet function at wavelet analysis?

Wavelets are functions which are form by two resulting coefficients. The detail coefficient and the scale coefficient. From wikipedia Wavelets are defined by the wavelet function ψ(t) (i.e. the mother wavelet) and scaling function φ(t) (also called father wavelet) in the time domain.

The wavelet function is in effect a band-pass filter and scaling it for each level halves its bandwidth.

The scaling and detail basically divide the signal into two applying a high-pass filter resulting into the detail coefficients - (which is the highest level of the transform) and a low-pass filter which results in the scaling coefficients - (which is the lowest level of the transform). This coefficients are calculated for each frequency of the discrete wavelet.

used to for example to extract the deterministic part of a signal and the stockastic part of the signal. They are also used through different procedures to denoise via thresholding the coefficients.


Gonzales&Woods is a good book.

since u need some other book for reference i may suggest

Author: S. Jayaraman, S. Esakkirajan And T. Veerakumar

Title of the Book: Digital Image Processing

Publisher: Tata McGraw - Hill Education Pvt. Ltd, New Delhi

Year: 2010

Edition: 3

This is also good book for beginners


getting different output when using FFTW3 Library than using Rosetta FFT

If the input is real valued than the output of the FFT is conjugate symmetric (Hermitian Symmetry). That means that the data in the second half of the output is redundant and the FFTW implementations takes advantage of that fact and only returns a complex vector of length $N/2+1$, i.e. 5 in your case. The Rosetta implementation just runs a standard complex FFT.
I suggest reading the documentation: http://www.fftw.org/fftw3_doc/One_002dDimensional-DFTs-of-Real-Data.html
Please note that your output is complex but you are only printing the real part.


Efficient calculation of correlation function every $N^{\rm th}$ point

If you have a signal of length $N$ you can calculate auto correlation function of length $K=2N-1$.You'll need perform $N$ multiplications for each value of correlation function $r_{xx}$ or in other words you need one correlator. So you need as many correlators as your $r_{xx}$ function length is. If you want to find correlation in the lag $0$ use one correlator. You can can calculate correlation of any lag $i \in (-N+1...N-1)$ via one correlator with signal properly shifted by $i$ times. From the other side $FFT$ method allows you to find $N$ taps of cyclic auto correlation function via only one $FFT$ operation. So it is not optimal if you need only a few correlation values.


Optimization of Filter Coefficients to Given Word Length

It is definitely a difficult mathematical problem (integer programming), but I don't see much research going on in the DSP community. I do not know about the maths community. There were some important papers in the early 1980's on the optimum design of FIR filters with quantized coefficients:

D. M. Kodek, "Design of Optimal Finite Wordlength FIR Digital Filters Using
Integer Programming Techniques," IEEE Trans. ASSP 28, 304-308 (1980).
D. M. Kodek and K. Steiglitz, "Comparison of Optimal and Local Search Methods
for Designing Finite Worldlength FIR Digital Filters," IEEE Trans. Circuits Systems 28, 28-32 (1981).

But as far as I know none of the published methods developed into a standard design method for filters with quantized coefficients, unlike for the design of filters with infinite precision coefficients (e.g., Parks McClellan).

I think there are two reasons for that: first, optimum methods use integer programming techniques, which are computationally inefficient, and which will only compute a local optimum that can be arbitrarily bad. These problems become aggravated for large filter lengths. Second, there is not so much need anymore for the optimum design of filters with quantized coefficients. This has to do with the fact that more and more systems are implemented in floating point, and - if fixed-point is needed - the word length is usually large enough such that suboptimal heuristic design methods are sufficient. Such methods iteratively design the filter assuming infinite precision, quantize the coefficients, check the filter properties and, if necessary, adapt the specs (including the filter length, error weighting, etc.), re-design the filter according to the new specifications, and so on, until the requirements are met.

For IIR filters, the optimum design for quantized filter coefficients is even more difficult than for FIR filters due to the rational transfer function. I haven't seen any useful optimization methods in the DSP literature, only heuristic methods such as the one described above. Note that different filter structures will result in different possible pole locations in the complex plane after coefficient quantization. The possible pole locations for the direct-form structure are relatively sparse at small pole angles, i.e., for low cut-off frequencies, whereas the pole locations of the coupled-form structure are evenly distributed on a rectangular grid.


If you merely want to shift spectral samples up in frequency, the process is fairly simple.  Let's say your DFT comprises the eight complex-valued samples X[0], X[1], X[2], X[3], X[4], X[5], X[6], X[7] and you want to shift them upward in frequency by two samples. The new shifted Y[m] spectral samples would be:

Y[0] = X[6]

Y[1] = X[7]

Y[2] = X[0]

Y[3] = X[1]

Y[4] = X[2]

Y[5] = X[3]

Y[6] = X[4]

Y[7] = X[5]

You'll notice that the shifting could be described as a "circular shift" because of the sample wrap-around behavior. If you want to shift your original X[m] sequence downward in frequency by three samples. The new shifted Y[m] spectral samples would be:

Y[0] = X[3]

Y[1] = X[4]

Y[2] = X[5]

Y[3] = X[6]

Y[4] = X[7]

Y[5] = X[0]

Y[6] = X[1]

Y[7] = X[2]


Why we need Laplacian for edge detection?

"Need" or "Want"?

Q1:  The second derivative avoids the problem of a gradually changing color being "greyish".

Q2:  Yes.

Here is one:

After advice about detecting focus quality of objects in a photo detected using YoloV3

It is based on how planar the color.  Very planar means no edge.

There are many, many other variations of these techniques possible.


Number of antenna to perform MUSIC algorithm

In the field of array signal processing, the signal model is
$$
\mathbf{x}(n) = [x_0(n), x_1(n), \ldots, x_{M-1}(n)]^T
$$
where $M$ denotes the number of sensor and $n$ is the discrete time.
While in the field of spectral estimation, the signal model is
$$
\mathbf{x}(n) = [x(n), x(n-1), \ldots, x(n-M+1)]^T
$$
where $M$ is the number of samples of the observed signal, and $x(n-m)$ represents the previous signal before $m$ samples.
The autocorrelation matrix is
$$
\mathbf{R}_{xx} = E\{ \mathbf{x}(n) \mathbf{x}^H(n) \}
$$
and it is usually replaced by its estimation $\mathbf{\hat{R}}_{xx}$ in practice.


You could check how "circular" a blob is, by computing the ratio between area and the square of perimeter.
For a circle that would be:

$A = \pi r^2 , P = 2 \pi r$,

which implies that

$ \frac{A}{P^2} = \frac {1}{4 \pi}  $

Probably those who have tail will have different ratio.
Both the area and the perimeter of the blob can be calculated in Matlab by using regionprops.

Another option [which could be interesting if you have a lot of data] would be computing many blob properties and training an SVM/Boosting/Other machine learning algorithm.


Fourier Transform Form: two sin components & a phase shift & a magnitude for only one term

Fourier Transform is a linear one, so you can make use of superposition principle:

$$ \mathscr{F} [\alpha x(t) + \beta y(t)] = \alpha \mathscr{F}[x(t)] + \beta \mathscr{F}[y(t)] $$

So for the first component $$x(t) = \sin \left( 2\pi \cdot 1000 \cdot t\right)$$

by definition:

$$\mathscr{F}\left[\sin(2\pi f_0 t + \phi) \right] = \dfrac{i}{2} \left[ e^{-i \phi}\delta(f+f_0) - e^{i \phi}\delta(f-f_0) \right] $$

you get:

$$ \mathscr{F}[x(t)]=\dfrac{i}{2} \left[ \delta(f+1000) - \delta(f-1000) \right]  $$

Second component is a sinusoid with shifted phase, so the complex exponent represents that:

$$y(t) = \dfrac{1}{2} \sin \left( 2\pi \cdot 2000 \cdot t + \dfrac{3\pi}{4} \right)$$

has following Fourier Transform:

$$\mathscr{F}[y(t)] = \dfrac{1}{2}\dfrac{i}{2} \left[ e^{\dfrac{-3\pi i}{4}}\delta(f+2000) - e^{\dfrac{3\pi i}{4}}\delta(f-2000) \right] $$

By summing both results you get the Fourier Transform of your signal.


Zero-pad before or after windowing for FFT

If done correctly, the order should not matter. You apply windowing to remove the jumps at the end of the signal segment by fading the values there to zero. Zero-padding then continues with zeros. Irregardless if padded or not, the windowing should only be applied to the original signal segment.

The wrong way would be to zero-pad and then apply windowing to the now longer signal over its full length, i.e., without regard to the position of the jumps. This would not remove the jumps and thus keep the background in the spectrum that is not related to the frequency content of the signal but the segmentation method.


what doest "hit or miss" in morphology do?

The Hit-or-Miss transform as the name suggest uses 2 structuring elements (SE) to identify structures which are specific to the foreground (first SE) and background (second SE).

Here, you have a good example for corner detection. The only difference is that the two SE are merged into 1: The 1 are part of the foreground SE in order to determine a specific shape in the pattern, and the 0 are part of the background SE. Consequently, the algorithm hits only if all the conditions on the 1 and 0 are respected, else it's a miss.

This principle is mainly used for skeleton computation, but it can be derived to build all the mathematical morphology basic operations.

Here are more explanations:
 - A complete video of a class.
 - Much more examples in this complete lesson.
 - StackOverflow