impulse response of a causal LTI system This is a difference equation to a causal LTI system:
$y[n] = ay[n - 1] + x[n] - a^Nx[n - N]$
Where N is a positive integer. I need to determine the impulse response of the system, so I have the equation:
$h[n] = ah[n-1] + \delta[n] - a^N\delta[n - N]$
h[n] is simple to find for the case where n < 0, n = 0, or n = 1. However, the response could be different for the case n = 2, depending on the value of N (in this case, N = 1 or N = 2):
$h[2] = ah[1] + 0 - a^1\delta[2-1] = a^2$
$h[2] = ah[1] + 0 - a^2\delta[2-2] = 0$
So which one is it? Intuition tells me the top one as the impulse response will be decaying, but we only know that N is positive.

HINT:
Just write down the values of the output signal $y[n]$ for $x[n]=\delta[n]$ for values of $n$ from $0$ to $N-1$ (you don't need any specific value for $N$, just use $n=0,1,2,\ldots$ and then you'll see what happens at $n=N-1$). Then figure out what happens at $n=N$, and what consequences this has on the output values for $n>N$. You may be surprised to find out that the given system has a finite length impulse response, even though it is recursive.


Non Equispaced / Non Uniform DFT Bandwidth I need to construct Fourier transform of non-equispaced data.

That is, I have signal $s(t)$, $t\in[0,T]$ sampled at non-equispaced points $t_k$, $k=0...N-1$ with sample values $s_k = s(t_k)$. For Fourier transform I use approximation of integral:
$$
   \hat S(\omega) = \int\limits_{-\infty}^\infty s(t)e^{-i\omega t}dt \approx \hat S_d(\omega) = \sum_{k=0}^{N-1}s_ke^{-i\omega t_k}\Delta t_k
\tag{1}
$$
where $\Delta t_k = (t_{k+1}-t_{k-1})/2$. As sampling points in frequency domain I choose $\omega_n = \frac{2\pi}{T}$.

My question is: as I can evaluate (1) for any $\omega$, what is the maximum $\omega$ such $\hat S_d(\omega)$ "adequately" represents $\hat S(\omega)$? What is the maximum of $n$ for $\omega_n$ I can use? For DFT we got Nyquist frequency. Do we have something similar for NDFT? Any references would be appreciated.

Note: I'm aware of such things as NDFT and NFFT. Though formula for NDFT, as presented in most papers, is
$$
\hat S_d(\omega_n) = \sum_{k=0}^{N-1}s_ke^{-i\omega_n t_k}
$$
I strongly believe that I need to use formula (1) as I'm trying to build a periodogram.

And I'm not interested in fast ways of computing NDFT yet, so I'm not considering NFFT.

Your formula isn't accurate.  Since you aren't trying for speed and internally  consider the interpolation between data points as a step/pulse then the formula should be.

$$
S(\omega)=\sum_{n=0}^{N-1}s_{k}sinc(\frac{\omega\triangle\left(t_{k}\right)}{2})e^{-wt_{k}}
$$
Having said that, what is "adequate"?  This is, more or less, subjective.  The original Nyquest criterion is based upon the result that any (freq < sam/2) (I think there is some small quibbling on that) can be perfectly reconstructed; if it and the sampling lasts long enough.  It says nothing about out of band signals; but the actual calculations will show aliasing.
It behooves you to generate the signals and interference that you envision and run it through your proposed algorithm. An alternative is to treat your system as a filter and try to "deconvolve" it. This is hazardous but can be done.
For a in depth look (with the hazards) read "Introduction to the Mathematics of Inversion in Remote Sensing and Indirect Measurements"  SIAM.
Now MRI/CAT scanners use various means of reconstructed signals.  All of the mathematicians seem to think that their technical success is an example of extreme luck; although I am sure that the designers might take exception to that.  You might try borrowing some of their mathematical techniques though.
Now in general instrumentation what we concentrate on is "noise temperature" and such.  This can be translated as: how far will your system signal/noise (or signal/distortion) be from the optimal number.  If you can get within 1db of the optimal S/N value then you should probably give up.  Of course this implies that you can calculate the optimal S/N.  In most instruments this can be done and this criterion gives you a stopping point; a point to give up and look for something else to do :)


Averaging data from 2 sensors Suppose we take a reading of a single rotating axis using 2 identical sensors. Intuitively, I think that averaging the outputs of both sensors will yield a less noisy result. But what if sensor A has a defect that makes it much noisier than sensor B. Will it still be wise to use the averaged data? Or will it be better to use only the data from sensor B?

If sensor A has a defect, the clear answer is to only use sensor B.
A preferred solution to minimize noise would be to do a weighted average based on the quality of each sensor, when that can be actively characterized. This can be easily done for the case the OP has presented of taking a reading from a single rotating axis.
The optimum combining for the general case of independent signals in the presence of white noise is to do a weighted average of each signal using the signal to noise ratio, which is referred to as maximal ratio combining. In practical application we determine the extent that we can use these assumptions, and our ability to extract an SNR metric depends also on how stationary the signals are. With consideration to this we arrive at a reasonable approach for the OP specifically.
A standard average of the result of two sensors given as x and y would be done as:
$$r_{avg} = \frac{x + y}{2} \tag{1} \label{1}$$
Alternatively, with maximal ratio combining, the two sensors could be combined as follows:
$$r_{opt} = \frac{\sigma_y^2 x + \sigma_x^2 y}{\sigma^2_x  + \sigma^2_y  } \tag{2} \label{2}$$
Where:
$\sigma^2_x$, $\sigma^2_y$: variance of each sensor's measurement
It is intuitively clear From \ref{2} how this would be superior, since we are now doing a weighted average with favoritism toward the less noisy sensor in the average of the two. [Addressing FAT32's good comment to an earlier version of this post; in addition to the variance estimate $\sigma^2$, a moving average estimating the mean $\mu$ of each sensor can be tracked for further robustness against erroneous measurements or a completely faulty sensor.]
Notice that when we have no further information as to the quality of each sensor, \ref{2} becomes \ref{1} when we otherwise then set $\sigma^2_x = \sigma^2_y$. To the extent that the measurement results are stationary, we are able to easily extract the first and second order statistics (mean and standard deviation as $\mu$ and $\sigma^2$) directly from the measurements of each sensor to then be able to apply maximal ratio combining. In many cases stationarity can be reasonably assumed over a finite time interval, allowing in the OP's case to easily estimate the mean and standard deviation from each sensor over the maximum time duration allowed by the dynamics of the experiment, bandwidth of sensors etc, and then combining as in \ref{2}.
For a very robust (but detailed) approach, see How to interpret Allan Deviation plot for gyroscope? for further details on using the two-sample variance to dynamically determine the best averaging time for the mean and variance estimates. Thus it can be envisioned where the measurement system is continuously monitoring a moving windowed mean, variance and two-sample variance of the sensors to provide the maximal ratio solution, where the window duration is updated based on minimizing the two sample variance. That level of complexity is likely unnecessary, as it may be sufficient to choose a moving window that is known to be a stationary interval and simply compute the mean and variance over that duration to get the result defined in \ref{2}; but shows an interesting robust approach of what else could be done.
In summary, Maximal Ratio Combining is superior solution over simple averaging. A simpler, suboptimum but effective strategy is to select the sensor with the highest SNR and ignore the other, or doing the simple average and disabling a sensor upon a threshold noise condition (which should be included regardless, such as a stuck sensor or excessive offset condition, to Fat32's points in the comments). Worst case would be to arbitrarily just average both for the very reason the OP brings up when it is entirely feasible that one of the sensors could become defective. A simple average would only make sense when there is high confidence both sensors are without defect and we believe to have similar performance.
This is similar to the strategies used in diversity receivers where two independent channels are optimally combined:
https://en.wikipedia.org/wiki/Diversity_combining
Also to include from the comments the link Fat32 found which refers to this approach as the "Fraser-Potter fixed-interval smoother", from the sensor community while I was exposed to equivalent "Maximal Ratio Combining" from the wireless communications community.
https://en.wikipedia.org/wiki/Sensor_fusion
As a demonstration per my friend Fat's request, below shows the results from a simulation of two sensor readings with a background trend, offset, and slope differential, where the statistics of the signals are not known apriori (everything is extracted from the signals themselves).
The top plot represents the raw measurements from the two sensors, with one sensor significantly noisy. The bottom plot shows the comparative results of simply averaging the reading from each sensor to produce an average result after each update (standard average) and the maximal ratio combining average where each of the two samples being averaged is first weighted by a variance estimate derived from the data from each sensor. The estimators used a moving window to measure the variance of each sensor over the past 500 samples.

An error signal was generated comparing the final results to the avg of the underlying sensor signals without the added noise and is plotted below showing the clear performance improvement in overall noise of the final result.

I'll also include the results of the Allan Deviation (root of the 2-sample variance) showing the significance behind it. I chose to use 500 samples from a review of the ADEV confirming for us in this particular case that statistical computations up to 500 samples will improve the overall estimates, anything beyond that or less than that will be degraded (this includes our estimates of the variance of each sensor for the maximal ratio combining to get each result as well as how long we can average those results over time for the final sensor result). This isn't a necessary step but something I do personally check in any statistical analysis and here demonstrates its clear utility. Basically this is saying that even though there is a background slope, given the noise level itself we can assume stationarity in this case over a 500 sample time interval. This means we would improve our noise estimates by using up to the past 500 samples for each measurement to determine the statistics for the weighting. (And using as much as 2000 samples would provide as good of an answer as 100 samples!). Also most significantly, note the clear difference between the standard average and the maximal ratio combining, in that the latter is better than the best sensor alone!
The horizontal axis for the plot below is in log samples, showing the minimum of the more dominant noise at a 500 sample averaging time.
The 30 second primer in how to read the ADEV plot: The horizontal axis indicates an averaging duration (here in samples) while the vertical axis is an indication of the standard deviation of the mean of a random variable after being averaged over that duration. This relates to our confidence in any estimator that may be used on the underlying data. For example, when determining the average value over $N$ samples, if the samples are random white noise and stationary, the standard deviation will go down by the rate of $1/\sqrt{N}$ which is exactly what we see in the plot for all curves out to 100 samples (and out to 500 samples for the blue curve of the noisiest signal). After that it bottoms out and turns up due to the underlying drift in the measurement. This tells us for all practical purposes we can assume the noise on the data is white and stationary out to the duration in time where the ADEV plot is still going down by the rate of $1/\sqrt{N}$, and that we would minimize the uncertainty (noise) of any estimators that include a mean as part of their computation (rms, variance, etc) by using data over that duration: very simply it tells us the true extent as to how long in time we can average for and keep improving a result; average over a longer duration and results actually get worst!   (In this case we could simply remove the linear drift to them take advantage of an even longer averaging time - but the minimum would typically occur due to the signal statistics being no longer stationary due to more random effects that cannot be predictably removed.)

TO DO: The Kalman Filter similarly determines optimum weights for the averaging between the two samples. Show how in this case the Kalman would result in the same weighting solution (it is the optimum combining for this case) and why, and show under what conditions specifically the use of the Kalman would be warranted...


Finding frequency respone of a differential/integral LTI system So suppose that we have an LTI system defined by the differential/integral equation below, where $x(t)$ and
$y(t)$ denote the system input and output, respectively. How would I find the frequency response of this systems?

I haven't been taught this yet and the text book has no examples on it, I'm really confused. I'd really appreciate the help. This is the system:

$$
\frac{d}{dt}y(t) + 2y(t) + \int_{-\infty}^t 3y(\tau)d\tau + 5\frac{d}{dt}x(t)-x(t) = 0
$$

Frequency response deals with the steady-state behavior. Differentiating with respect to $t$ yields:

$$\frac{d^2y}{dt^2}+2\frac{dy}{dt}+3y(t)=\frac{dx}{dt}-5\frac{d^2x}{dt^2}$$

The Laplace transform is $$(s^2+2s+3)Y(s)=(s-5s^2)X(s)$$
and we have

$$H(s)=\frac{Y(s)}{X(s)}=\frac{s-5s^2}{s^2+2s+3}$$
I think you can easily do the rest.


Methods to quantify randomness (or complexity) in a signal What are the methods through which we can quantify the randomness or complexity in a given signal. I know spectral flatness measure (geometric to arithmetic means) is one way to do it, but what are other common techniques?

There are two different concepts:

If you think as your signal as a single random variable $X$ that is emitting values, then what you want is to calculate the Entropy of the random variable

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

If you are considering the entire random signal or stochastic process, then you have to estimate the autocorrelation function. The most random signal possible is White noise, in which every sample of the signal is not correlated with any other except itself.

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


How Do I Plot Color Distribution of HSV Image in MATLAB I have a dataset of images representing human skin. How do I plot the skin color distribution in HSV color space and YCBCR color space. However, I'd like ignore the "value" or "brightness" in both color spaces and analyze rest of the layers.

edit

The below is the code, I found to plot the distribution for ycbcr but don't understand what it means. Can anyone explain it please?

chroma = zeros(256);
cb = imycc(:,:,2);
cb = reshape(cb, 1, numel(cb));
cb = round(cb);

cr = imycc(:,:,3);
cr = reshape(cr, 1, numel(cr));
cr = round(cr);

for i = 1:length(cb)
   chroma(cb(i), cr(i)) = chroma(cb(i), cr(i)) + 1;
end
surf(chroma)

Have a look at File Exchange Submission - RGB/HSV Distribution on  visualization of RGB/HSV distribution in a given image:

hsv_distribution(inputImage, 5)

It should do what you're after as the result is like:


Coefficients for half-order filter var white = Math.random() * 2 - 1;

// Filter to approximate 10 dB per decade
b0 = 0.99886 * b0 + white * 0.0555179;
b1 = 0.99332 * b1 + white * 0.0750759;
b2 = 0.96900 * b2 + white * 0.1538520;
b3 = 0.86650 * b3 + white * 0.3104856;
b4 = 0.55000 * b4 + white * 0.5329522;
b5 = -0.76162 * b5 - white * 0.0168980;
var sample = (b0 + b1 + b2 + b3 + b4 + b5 + b6 + white * 0.5362) * 0.15; // Roughly compensate for gain
b6 = white * 0.115926;

Above is some code that approximates a 10 dB/decade filter, where "white" represents a random number between -1 and 1, producing white noise, while "sample" is the actual value that is passed as a sample to the output buffer.

My question is, how are the coefficients 0.99886, 0.0555179, et cetera actually determined? Including the 0.5362 and 0.115926

Allan Herriman notes on the music-dsp mailing list that the filter, created by Paul Kellett and named the "instrumentation grade" pinking filter for white noise, is not optimal in (my words) minimum absolute error sense in the magnitude frequency response. I don't know if it should be, perhaps it would be more desirable for it to be that on dB scale. My guess is that the coefficients are handcrafted. Also I don't know if the filter tries to approximate a "half-order" filter or to have an approximately 10 dB / decade slope over the audio band for pink noise generation from white noise, probably the latter.


How to avoid overtones in audio signal I'm trying to make a program that catch and audio signal of a bass and show the notes that have been played. My program:

Gets the sound and put it in a double data type array. (48 000Hz, 16bits) (1 frame in array per sample)
Splits the array into many arrays, with 4~5 arrays in a second approx.
Applies FFT algorithm to every array and left the same number of array's that have been created in the last step with frequency domain instead of time domain.
Obtains the highest number of each array. That's the note that I have played with the bass.

This works and get the frequency when I don't split the array, for most of the times, but when I split it doesn't work even for the half of the results, maybe I have to use another way for get several notes instead of splitting. But the results are good too, so I think the problem is that I have to apply some filters.

The first one has to be a filter for the range of the bass (41 - 784), assuming that I'm not playing anything that time if it's out of that range.

Many times, if I play 41Hz (E1), the program interpret that I'm playing 82Hz (E2), or even a higher harmonic, and I don't know what to do with that. I could turn into 0Hz every note that has been played after the same note or one of its harmonics, but then how could I play the same note several times.

Some advice for deal with this problems?
Should I split the original array into more/less array's for get more accuracy?
Which filters should I apply on my frequency or time domain array for get all the notes that have been played, with the right frequency and avoiding their harmonics?

A good stringed instrument music pitch detection/estimation algorithm will do the opposite, e.g. it will not ignore overtones.  Instead it will pay attention to the harmonics, specifically the harmonic train and its spacing, as this is a stronger indication of human perceived pitch than spectral content at the fundamental frequency (which could be nearly or completely missing, depending on the microphone and channel characteristics).

If you want to work in the frequency domain, than look into cepstrum/cepstral analysis and/or the Harmonic Product Spectrum algorithm, which can estimate pitch from spectral periodicity.  Or you can instead work in the time domain using weighted/interpolated autocorrelation, ASDF, AMDF, etc. type algorithms for pitch estimation.

If you slice your time windows too small, then there will be very few full cycles at he fundamental pitch frequency, which causes FFT frequency estimation to have very poor resolution.  A parametric estimator within and across windows may work better than a bare FFT.


Why does changing the frequency of a single component change total sum of amplitudes in my FFT? I have the following code which generates a single simple sine wave with amplitude $1$ and then calculates the Fourier transform, normalises it, and finally takes the sum of all the components.

However I've observed something strange, if I generate a $10\textrm{ Hz}$ wave and a $100\textrm{ Hz}$ wave. The amplitudes are pretty damn similar as expected and they are  0.9999483519908307 and 0.9997854546246567 respectively.

However the total sum of all the components after normalisation are 1.01723486798 and 1.17263231047 respectively, which is approximately a $15\%$ difference! Why is this occurring?

Below are the plots of the spectrums and the python code below.

mport matplotlib.pyplot as plt
import numpy as np
from decimal import Decimal

#number of samples
iNumSamples=10000

#period of samples (i.e. sampling at 10Khz)
Ts=1/iNumSamples

#time axis used to sample 10,000 times for a second
t = np.linspace(0, 1, iNumSamples)

#10Hz Signal
signal1=np.sin(t*20*np.pi)

#100Hz Signal
signal2=np.sin(t*200*np.pi)

#creates a new figure object
freqFig=plt.figure(2)
#change title
freqFig.canvas.set_window_title('Frequency domain')

#fast fourier transform of 10Hz signal
absoluteValues1=abs(np.fft.rfft(signal1))

#fast fourier transform of 100Hz
absoluteValues2=abs(np.fft.rfft(signal2))

#find the frequencies of the rtansform
frequencies=np.fft.rfftfreq(iNumSamples, d=Ts)

##Normalises (rescales) the fft in a sensible way so that the corresponding freq component of a sine wave is 1 if the sine wave amplitude is 1
##non zero freq amplitudes (even when using rfft) are halfed as there are negative and positive components, hence first multiple by 2

#first create empty array of correct size
normalisedValues1=np.zeros(absoluteValues1.size, dtype=Decimal)
normalisedValues2=np.zeros(absoluteValues2.size, dtype=Decimal)

#normalise 0 component
normalisedValues1[0]=absoluteValues1[0]/iNumSamples
normalisedValues2[0]=absoluteValues2[0]/iNumSamples
#normalise non zero components
normalisedValues1[1:]=absoluteValues1[1:]*(2/iNumSamples)
normalisedValues2[1:]=absoluteValues2[1:]*(2/iNumSamples)

print("peak value signal 1:")
print(np.amax(normalisedValues1))
print("peak value signal 2:")
print(np.amax(normalisedValues2))

print("total sum of all componenets in spectrum of 10 Hz")
print(np.sum(normalisedValues1))

print("total sum of all componenets in spectrum of 100 Hz")
print(np.sum(normalisedValues2))

#creates a new figure object
freqFig=plt.figure(1)
#change title
freqFig.canvas.set_window_title('Frequency domain')
plt.plot(frequencies,normalisedValues1)

#creates a new figure object
freqFig=plt.figure(2)
#change title
freqFig.canvas.set_window_title('Frequency domain')
plt.plot(frequencies,normalisedValues2)

plt.show()

That's a subtle sampling issues and you will get your desired result by replacing this line:

t = np.linspace(0, 1, iNumSamples)

with the following:

t = np.linspace(0, 1-Ts, iNumSamples)

That being said, it's really a virtual sucess and you would never rely on such perfectly matching numbers in practice...


Wiener filter optimal coefficients and the input signal I have this question:
"Show that for the optimal filter coefficients the error $e[n]$ is orthogonal to the input signal $w[n]$."
anybody have an idea?

For conveniece I'm going to drop the time indexes. Let $e=d-\mathbf{a}^T\mathbf{x}$, where $d$ is the desired signal and $\mathbf{a}$ weights the elements of the observations/input vector $\mathbf{x}$. The goal is to minimize expected squared error with respect to the vector $\mathbf{a}$.

$$\min J=\min E[e^2],$$

where $E$ is the expectation operator.
To find the minimum we take the derivative of the objective function w.r.t. $\mathbf{a}$ and set it equal to zero. Therefore
$$
\mathbf{g}=\frac{\partial E[e^2]}{\partial \mathbf{a}}=E[2e\frac{\partial e}{\partial \mathbf{a}}]=0,$$

but from the defintion of $e$ earlier, we have
$$ \frac{\partial e}{\partial \mathbf{a}}= -\mathbf{x}$$
Therefore
$$ E[2e\frac{\partial e}{\partial \mathbf{a}}]= -E[2e\mathbf{x}]=0. $$ The factor 2 can be removed and thus $E[e\mathbf{x}]=0$, which means the error $e$ is orthogonal to the input signal $\mathbf{x}$.


What does correcting IQ do? Some software (for example SDR Sharp) has a "Correct IQ" option.  This is used to get rid of the DC spike in the middle of the display.  What does this setting do to the signal?  What is wrong with the IQ values that needs to be corrected?

Here are some graphics I have to add to the other good answers, as well as typical correction schemes for DC Offset and IQ Imbalance.
DC OFFSET
Below is an example DC Offset on a 16 QAM signal (with some added AWGN). In this example the offset if shifted by an amount such that the power in the fixed vector representing the shift is 10 dB below the total power in the signal. Notice that the shift can be at any magnitude and phase and the fact that it does not change with time makes it "DC". Also notice that even though it is 10 dB BELOW the total power in the signal, it appears significantly higher in the spectrum plot which may confuse some observers: this is because the waveform's power is spread over the modulation bandwidth, so the power in each FFT bin is lower than the power in the DC offset.

Below is a simple approach for eliminating DC offset. This filter is also introduced in this post: Transfer function of second order notch filter, where the gain constant $\alpha$ controls the bandwidth, getting tighter as $\alpha$ approaches 1.

IQ IMBALANCE
IQ Imbalance is demonstrated on the same 16 QAM signal in the next figure below, showing how it can result from either amplitude imbalance or phase imbalance or both. With amplitude imbalance, the scale of one axis is changed from the other, with phase imbalance, the axis are no longer in perfect quadrature.

The correction for IQ Imbalance is done as in the figure below with an $\alpha$ gain factor shown for amplitude imbalance and a $\beta$ gain factor shown for quadrature imbalance.


Difference between these two digital down conversion methods I can think of two possible methods of digital down conversion. One seems superior to the the other for most cases, but I'd like to get some DSP experts' comments on the practical differences between these two.

For both methods described below, I start with the assumption that the $x(t)$ input is already constrained to the band of interest, to avoid the consideration of antialiasing. Antialiasing filters may be necessary in practical applications.

In method A, the in phase component is generated by mixing $x(t)$ with $\cos(\omega t)$ and low pass filtering out the undesired high frequency images. Since there is excess bandwidth in the signal you can use a decimating filter instead of a straight lowpass, shown in my figure below as a blue box around a low pass and a decimator. The quadrature component is similarly obtained by mixing with $\sin(\omega t)$ and then filtering. As a formality I show multiplying the quadrature component by $j$ and adding it to the in phase component.

In method B, convolution is used to perform the discrete Hilbert transform on $x(t)$. This is the reconstruction of the missing quadrature component in $x(t)$. Downconversion is accomplished by mixing the resultant complex signal with a complex sinusoid, no unwanted images generated! Depending on the bandwidth of the baseband signal, it may be smart to add a decimating filter to reduce the data rate, but the output here is usable without that filter. That is, it is not aliased (provided the input signal is limited to the appropriate band).

The primary constraint that I can see in method A is that you can't downconvert by small amounts relative to your bandwidth. Say, for example, you had a 1kHz bandwidth in your signal of interest. If you want to downconvert by 100Hz, method A will not allow for that. The band of the images generated in mixing will intersect the band of the signal of interest. I've seen vector analyzers with analog front ends that can use double down conversion to reject images. A similar method could be used in DSP of course. The problem isn't unsolvable, it just isn't solved by method A.

Method B appears to work around the problems with method A. The construction of a complex signal from the real signal essentially eliminates all negative frequency components, and mixing with a complex sinusoid ensures no images are created, so the only signal present in the output is the signal desired. No filtering is required. However, the convolution required to do the discrete Hilbert transform could be expensive. I would generally choose method A.

My question is: are there advantages / drawbacks of these methods I haven't thought of? What are some practical considerations I should keep in mind before choosing a method?

As you already know, the two methods are theoretically identical if there's no noise. If there is out-of-band noise then the low pass filters in method A will further suppress the out-of-band noise, whereas there's no such noise suppression with method B. This would be one advantage of method A over method B.

Note, however, that the phase splitter in method B (which is the complex-valued filter with a wire in the real part and a Hilbert transformer in the imaginary part) could be replaced by a complex valued band pass filter filtering out the positive in-band frequencies. With such a filter, both methods would be equivalent, even with out-of-band noise.

I do not think that the need to implement a Hilbert transformer is a disadvantage of method B, because the Hilbert transform does not need to be approximated over the whole band, just over the signal bandwidth.


Bilinear interpolation implemented by convolution I read the paper Deep Feature Flow for Video Recognition https://arxiv.org/abs/1611.07715.

In Sec.3, the author implements bilinear interpolation like this:

$$f_i^c(p)=\sum\limits_{q}G(q,p+\delta p)f_k^c(q) \tag{1}$$

Where $q$ is the point from the source image, and $p$ is the points on the target image. $\delta p$ is the distance the point moved each point $p$ (not $\delta \bullet p$). $G$ is defined as

$$G(q,p+\delta p)=g(q_x,p_x+\delta p_x)g(q_y,p_y+\delta p_y)\tag{2}$$

And the bilinear interpolation is defined in wiki as:

$$f(x,y)\approx {\frac {y_{2}-y}{y_{2}-y_{1}}}f(x,y_{1})+{\frac {y-y_{1}}{y_{2}-y_{1}}}f(x,y_{2})\tag{3}$$

I think the operation $(1)$ and $(3)$ is equivalent. How can I derive the filter $(1)$ from $(3)$?

Honestly I haven't read that article you linked to, but as long as you want a convolution kernel for 2D bilinear interpolation, then the following should help.

Bilinear interpolation gives a crude result which can be sufficient in case the application does not require a perfect output otherwise. Then the following simple triangular 2D convolution kernel implements a bilinear interpolation of the ratio (after compression by $M$) given by $r = \frac{L}{M} > 1$.

L = 9;          % upsampling ratio
M = 5;          % downsampling ratio (M < L assumed)

h   = (1/L)*conv(ones(1,L),ones(1,L));  % 1D linear interpolator.
hBL = h'*h;                             % 2D bi-linear interpolator.

I = im2double(imread('Cameraman.tif')); % Read some image into double
S = size(I);                            % size of image

Ie = zeros( L*S(1), L*S(2) ) ;          % expand the image by stuffing zeros
Ie( 1:L:end , 1:L:end ) = I;            % assign the original pixels

Iup = conv2(Ie,hBL);                    % linear upsample (interpolate) by L
Iint = Iup(1:M:end,1:M:end);            % downsample by M

figure,imshow(I);title('original image');
figure,imshow(Iint);title('upsampled by L/M times');


How to get histogram data object from matlab Lets say I have a matrix x=[ 1 2 1 2 1 2 1 2 3 4 5 ].
To look at its histogram, I can do h=hist(x).

Now, h with retrieve a matrix consisting only the number of occurrences and does not store the original value to which it occurred.

What I want is something like a function which takes a value from x and returns number of occurrences of it. Having said that, what one thing histeq does should we admire is, it automatically scales nearest values according!

How should solve this issue? How exactly people do it?

My reason of interest is in images:

Lets say I have an image. I want to find all number of occurrences of a chrominance value of image.

%Let's say:
X = [1 2 3 4 4 4 3 2 1 5 10];
Num = input("Please enter the number");
temp = X == Num;

%let's assume Num in this case is 4
%temp = [0 0 0 1 1 1 0 0 0 0 0]

ans = sum(temp);

ans = 3

This is a quick and an effective way to find the number of occurrences of a particular value.


Detect almost circular shapes in RGB image I want to extract circular trees from this rgb image, I've tried circular hough and segmentation techniques but in vain. Any help?

The image itself is somehow blurred a little. Whole edges of trees are not visible. The only clue of tree existence seems like shadows. So I suggest using shadow information and watershed transform and local thresholding. I have done similar exercise in here.  The result is not that perfect but it will give some idea. Here are the steps. 1. Find shadows, 2. Expand shadows towards sun direction 3. Create shadow mask 4. Blur the image 5. Watershed transform 6. Threshold each watershed segment individually 7. Fit circle with Hough transform

im=imread('P0Yvl.png');
imOrig=im;

% Grayscale
im=rgb2gray(im);

% Shadows
shadow=im<60;
[xs,ys]=find(shadow);

% Grow shadow towards sun direction

% Create structural element for dilation
str_element_initial=flip(eye(25));
str_element_initial(:,1:12)=0;

shadow2=zeros(size(str_element_initial));
for i=1:7
    shadow2 = shadow2 | imresize(imrotate(str_element_initial,40-i*10),size(str_element_initial))>.1;
end

mask_shadow=imdilate(shadow,shadow2);
mask_shadow=mask_shadow-shadow;

% Plot shadow mask and initial shadow locations
figure,imagesc(mask_shadow);axis image;
hold on, plot(ys,xs,'r.')

% Blur image
im2=imfilter(im,fspecial('gaussian',[21 21],3));

% Watershed transform
L = watershed(max(im2(:))-im2);
[xw,yw]=find(L==0);

figure,imagesc(imOrig),axis image
hold on, plot(yw,xw,'r.','MarkerSize',1)

% Threshold each segment individually
tmp=zeros(size(im));

for i=1:max(L(:))
    mask=L==i;
    ind=find(mask);
    thr =multithresh(im(ind),1);
    tmp(ind)=im(ind)>thr*.8;
end

% Remove detections that don't have shadow
trees=tmp & mask_shadow;
trees=bwareaopen(trees,10);

% Find circles with Hough transform
[centers, radii, metric] = imfindcircles(trees,[5 25]);

% Result
figure,imagesc(imOrig),axis image
viscircles(centers, radii,'EnhanceVisibility',0,'LineStyle','-','LineWidth',1);


what is matlab function is taking as sign of second order IIR filter? I am getting confused from the implementation regarding matlab function filter([b0,b1,b2],[a1,a2],x) and other function similar as freqz. I do not understand which sign I must put in a coefficients, since there is two different way of writing the transfer function for the same second order filter as you see on the these pic below

when for example the besself function generate the coefficients, does it generate them as I write them $H(z)$ or I must write $1-a_0z^{-1}-a_2z^{-2}$?

The filter function assumes this form (for your case):

$$
\tilde{H}(z) = \frac{b_0+b_1z^{-1}+b_2z^{-2}}{1+a_1z^{-1}+a_2z^{-2}}
$$

where the $1$ on the bottom is from the assumption that $a_0=1$, and call it $\tilde{H}(z)$ so not to be confused with your $H(z)$. Since you have written the denominator with the signs of $a_1$ and $a_2$ inverted, then you just have to fix that before using filter.

Suppose someone wrote the transfer function from your question down $\big(H(z)\big)$, handed you the coefficients, and said use MATLAB to find the output. The only thing you need to do before is to invert $a_1$ and $a_2$:
filter(..., [a0, -a1, -a2], ...).

When using besself, it returns the coefficients in the form according to the transfer function $\tilde{H}(z)$, and you can then directly use those coefficients.


How to simplify this sequence So I had this question (f) at my last exam that I failed and I'm now studying for my re-exam, the problem is that I still don't know what the answer is supposed to be in this question (f):

Do you have any suggestions how to solve it?

HINT:

The only thing you need to know is that any sequence convolved with a (shifted) delta impulse is a shifted version of itself:

$$f[n]*\delta[n-k]=f[n-k]$$

Apply this to the expression $u[n]*(\delta[n]-\delta[n-1])$, and realize that what you're left with is a very simple sequence. (The result of (a) will be helpful.)

Then convolve this very simple sequence with $x[n]$ to get the final result.


Is there a way to derive an FIR filter using an IIR filter? I know there's a thread similar to this one here, but the OP is asking the reverse of what I'm trying to find here. I've done some research on the web with very few sources coming up with actual solutions to this problem. What techniques are used to give an approximation of an FIR filter given one or more IIR filters with say the same order?

Approximating the frequency response of an IIR filter or physical process using an FIR filter is useful in learning control. It is quite common to do FIR filter design based on frequency response specifications. You probably want to check out two standard papers on the subject:

[1] J. H. McClellan, T. W. Parks, and L. R. Rabiner, “A computer program for designing optimum FIR linear phase digital filters,” IEEE Trans. Audio Electroacoust., vol. 21, no. 6, pp. 506–526, 1973.

[2] L. R. Rabiner, “Techniques for Designing Finite-Duration Impulse-Response Digital Filters,” IEEE Trans. Commun. Technol., vol. 19, no. 2, pp. 188–195, Apr. 1971.

Broadly, you either do windowed direct sampling of your desired frequency response, or you use one of several optimization methods to achieve similar results. If you disregard the linear phase delay in an FIR, you can practically make the IIR and FIR responses identical, if the FIR filter order is high enough.

As an elaboration on one of the other answers given; if you have an IIR filter $G(z^{-1})$, then you can do FIR filter design by frequency sampling by taking $M$ samples of the frequency response of $G(z^{-1})$, denoted $\widehat{G}(k)$, and then taking the inverse discrete Fourier transform (IDFT) of $\widehat{G}(k)$. The unit impulse response $g(n)$ of $\widehat{G}(k)$ is
\begin{align*}
 g(n) = \frac{1}{M} \sum\limits_{k=0}^{M-1} \widehat{G}(k) \text{e}^{j \frac{2 \pi k n}{M}} \; ,
\end{align*}
where $n \in [0, M-1] \cap \mathbb{N}_{0}$. The FIR filter is then expressed in the $z$-domain as
\begin{multline*}
 F(z^{-1}) = g(0) + g(1) z^{-1} + ... + g(M-1) z^{-M+1}
 = \sum_{n=0}^{M-1} g(n) z^{-n} \; .
\end{multline*}

The frequency-sampling method results in a unit impulse response which has been convoluted with a rectangular window of the same length in the frequency domain. The frequency response of $F(z^{-1})$ is therefore affected by the large side-lobes of the rectangular window. As a result, the approximation error of $F(z^{-1})$ is large between the frequency samples. This can be alleviated by the use of a window that do not contain abrupt discontinuities in the time domain, and thus have small side-lobes in the frequency domain, i.e., the window smooths the frequency response of $F(z^{-1})$.

A windowed FIR filter $\tilde{h}(n)$ is created from an un-windowed FIR filter $h(n)$ as
\begin{align*}
 \tilde{h}(n) = w(n) h(n)
\end{align*}
where $w(n)$ is a window function which is non-zero only for $n \in [0, M-1] \cap \mathbb{N}_{0}$. The frequency-domain representation of the window function $W(k)$ is found as
\begin{multline*}
 W(k) = \sum\limits_{n=0}^{M-1} w(n-M/2) \textrm{e}^{-j \frac{2 \pi k n}{M}}
 = \left[ \sum\limits_{n=0}^{M-1} w(n) \textrm{e}^{-j \frac{2 \pi k n}{M}} \right] \textrm{e}^{-j \frac{2 \pi k}{M} \frac{M}{2} } \; ,
\end{multline*}
where the term $\textrm{e}^{-j (2 \pi k / M) (M/2) }$ comes from the fact that the rectangular window is not centered around $n=0$, but is time-shifted to be centered around $n=M/2$. This phase term will cause distortion of $h(n)$, unless $h(n)$ is also phase-shifted to compensate. The unit impulse response $g(n)$ is therefore phase-shifted before windowing. Due to the circular shift property of the DFT, this can be done by rearranging $g(n)$ such that
\begin{equation*}
 \bar{g}\left( n \right) =
 \begin{cases}
  g\left( n + M/2 \right) , & \hspace{-0.6em} n = 0,1, ..., \frac{M}{2} - 1 \\
  g\left( n - M/2 \right) , & \hspace{-0.6em} n = \frac{M}{2},\frac{M}{2}+1, ..., M-1
 \end{cases}
\end{equation*}
for the case when $M$ is even. The response is then represented by the FIR filter
\begin{equation*}
 \bar{F}(z^{-1}) = \sum_{n=0}^{M-1} \bar{g}(n) z^{-n} = z^{-M/2} F(z^{-1})
\end{equation*}
which is $F(z^{-1})$ delayed by $M/2$ steps. Applying the window $w(n)$ to the time-shifted impulse response $\bar{g}(n)$,
\begin{equation*}
 \tilde{g}(n) = w(n) \bar{g}(n) \; ,
\end{equation*}
the filter
\begin{equation*}
 \tilde{F}(z^{-1}) = W(z^{-1})*\left[ z^{-M/2} F(z^{-1}) \right]
\end{equation*}
is obtained. Now, $G^{-1}(z^{-1}) \left[ W(z^{-1})*F(z^{-1}) \right] \approx 1$ if the FIR filter is accurate. Note that the phase due to $z^{-M/2}$ is taken out.


Combined signal folding and shifting Let $g[n] = f[-n]$(n is the time index). Then what will be expression for $g[n+k]$.

Is $g[n+k] = f[-n + k]$ or $g[n+k] = f[-n - k]$. Which one is true and Why?

The statement $g[n]=f[-n]$ is meaningless unless you assign some numerical value to $n$, e.g. $g[3]=f[-3]$, and if you do say that $n=3$, then the equality
$g[3] = f[-3]$ tells you nothing about the relationship between $g[n+k]=g[3+k]$
to the value taken on $f[\cdot]$ for any choice of argument.

What you really meant to say is that

  $$\text{For all integers } n, g[n]=f[-n]$$

which asserts that the equality holds for all choices of $n$, and in
this case, it follows by definition that $g[n+k] = f[-n-k]$;
you just need to negate the argument $n+k$ of $g$ (to get $-n-k$)
and use it (it meaning $-n-k$) as the argument of $f$ exactly as
the definition says you should.  The
complicated arguments in another answer that claim that $g[n+k]=f[-n+k]$
are wrong, as discussed in the comments following that answer.


Checking Parseval's Theorem for Gaussian Signal by Using Scipy I'm trying to check Parseval's theorm for Gaussian signal. It's well known that fourier transform of $\exp(-t^2)$ is $\sqrt{\pi}\exp(-\pi^2 k^2)$. So I implement it by using quad and simps. I think that the one is continuous integral, the other is integration which uses sample.

Code

from scipy import integrate
import numpy as np

def f(x):
    return (np.exp(-x**2))**2
def F(k):
    return ((np.pi**0.5)*np.exp((-np.pi**2)*(k**2)))**2

a=integrate.quad(f,-np.inf,np.inf)
print(a)
#>>>(1.2533141373155017, 4.4674503165883495e-09)
b=integrate.quad(F,-np.inf,np.inf)
print(b)
#>>>(1.2533141373155066, 1.592743224014847e-08)

a=b so above code follows Parseval's Theorm. The problem is on following code.

import scipy.fftpack as fft

N=50000
t  = np.linspace(-1000,1000, N)
h=np.exp(-t**2)
H=2*np.abs(fft.fftshift(fft.fft(h)/N))
freq=fft.fftshift(fft.fftfreq(H.shape[0],t[1]-t[0]))

S_h=integrate.simps(h**2,t)
print(S_h)
>>>1.25331413732

S_H=integrate.simps(H**2,freq)
print(S_H)
>>>1.25326400525e-06

S_h is not equals to S_H. What is the problem? How can I fix it??

There is a simple idea behind Parseval's Theorem.
Unitary transform preserves the $ {L}_{2} $ Norm.
The Fourier Transform (All its variants) are indeed Unitary when the correct factor is in place.

The problem is with the way you sum the elements and factor them.
What you need to factor by $ \frac{1}{N} $ is the sum of squared absolute value of the series.

Have a look on the following MATLAB Code:

numSamples = 50000;
vT = linspace(-1000, 1000, numSamples);
vH = exp(-(vT .* vT));

vK = fft(vH);

% Summation by Parseval Theorem

sum(vH .* vH)
sum(vK .* conj(vK)) / numSamples

In MATLAB it yields

ans =

   31.3322

ans =

   31.3322

Why is that?
Pay attention that the fft() applied by MATLAB and Python is not the Unitary form.
For the unitary form it should be factored by $ \frac{1}{\sqrt{N}} $.

Yet in Parseval we see:

$$ \sum_{n = 0}^{N - 1} {x}^{2} \left[ n \right] = \frac{1}{N} \sum_{k = 0}^{N - 1} {\left| X \left[ k \right] \right|}^{2} $$

The factor $ \frac{1}{N} $ is exactly $ {\left( \frac{1}{\sqrt{N}} \right)}^{2}$.


Working with a sound's magnitude instead of amplitude I'm working on a project, where we're recording sound with a piezo-disc which looks a little something like this:

Now, unless we're doing something horribly horribly wrong, I've discovered that we're only recording the magnitude of sound, since a piezo-disc emits a response to pressure, regardless of direction.
This is a sound-recognition project and previously, we'd been working with smartphone microphones. Back then, we could MFCC to perform feature extraction, yet now it looks like some things will break:

I fear that the pre-emphasis step in particular might lead to faulty results.
What kind of an approach could we follow?

Is there a way to "guesstimate" the amplitude from the magnitude?
Assuming that I do skip pre-emphasis, I don't expect the Fourier transform of the magnitude over time to differ from that of the Fourier transform of the amplitude over time. But am I better of skipping it, or fine-tuning it?
I should just ignore all of these concerns and just continue to take MFCC's as usual.

Any other ideas?

There is no reason why your piezo shouldn't be able to produce a bipolar output, if you use proper biasing and/or preamp. See for example https://www.homemade-circuits.com/diy-contact-mic-circuit/
Working with the magnitude only is probably a non-starter. $y = |x|$ is a highly non-linear operation and will dramatically change your spectrum. For example a sine wave at $f_0$ has just a single spectral line. The spectrum of the magnitude has a strong DC bias, lots of even harmonics $2f_0, 4f_0 ...$ but no content at $f_0$ at all anymore


Continuous-time vs Discrete-time Fourier transform case 1) to calculate the Fourier transform of discrete-time signal(sampled signal) we use Discrete-time Fourier transform.

but my question is:

case 2) if I consider that discrete-time signal as continuous-time signal which is multiplied by Dirac comb(impulse train) then we produced a sampled signal but in continuous mode
now we apply continuous Fourier transform to this signal.

notice that in both case we have the same sampled signal(the only difference: in first case signal is in discrete form and in second signal is in continuous form)

the result of "case 1" is equal to result of "case 2"?
if not,what does cause it?

This is pretty straightforward to show. Let $x[n]$ be a sampled version of a continuous-time signal $x_c(t)$:

$$x[n]=x_c(nT)\tag{1}$$

The DTFT of $x[n]$ is defined by

$$X(e^{j\omega})=\sum_{n=-\infty}^{\infty}x[n]e^{-jn\omega}\tag{2}$$

The CTFT of the signal

$$x_c(t)\cdot\sum_{n=-\infty}^{\infty}\delta(t-nT)=\sum_{n=-\infty}^{\infty}x_c(nT)\delta(t-nT)\tag{3}$$

is given by

$$\begin{align}\mathscr{F}\left\{\sum_{n=-\infty}^{\infty}x_c(nT)\delta(t-nT)\right\}&=\sum_{n=-\infty}^{\infty}x_c(nT)\mathscr{F}\left\{\delta(t-nT)\right\}\\&=\sum_{n=-\infty}^{\infty}x_c(nT)e^{-jn\Omega T}\tag{4}\end{align}$$

which equals $(2)$ with $x[n]=x_c(nT)$ and $\omega=\Omega T$ being the normalized angular frequency.

And now I found this question, which is pretty similar but also asks about some background.


Decoding final bits with viterbi I'm writing a viterbi decoder and I'm rather confused about the correct process of decoding the final k bits for a kth order convolutional filter.

For the bulk of the message, my process considers the 0-successor and 1-successor of every state, that is, the resulting shift register state found by shifting a 0 in or a 1 in, respectively. Both of these successors share a least-error ancestor found by the dynamic programming step in the Viterbi algorithm.

For the last k bits, we know that the encoder shifted in 0s in order to flush out the payload and reset to the 0 state. I can confirm this is the case because I've also written the encoder. So, when I get to this part of the message, my decoder only considers 0-successors. I think this is all I can do with this information, but I feel like I'm possibly missing some extra info that would help recover the message.

What I'm discovering in reality is that unrecoverable errors are far, far more likely to occur in the final k bits of the message leading up to these 0 flush bits. For example, if my decoded message is 256 bytes long, and I inject errors at a rate of 3%, I find that 15% of the unrecovered errors occur in the last byte.

To some extent, I might expect this to be the case as there aren't additional bits left to help the Viterbi best path converge. I'm using a convergence length of 5k, or 35 decoded bits for a 7th order filter. Should I just expect that without a padding of 35 bits, the end of my message is more likely to receive errors? Or should I keep looking for a bug in my decoder implementation? In practice, is it common to add additional payload bits to help viterbi converge?

You do not need any additional bits to help the decoder converge.  There should be no need for convergence at the boundaries, because the start and end states are completely known (the all zeros state) so there is no guessing involved.

If you are getting more errors at the end it is likely due to a bug in the decoder.


Peak location estimation for a triangle wave in AGWN? Is there a more optimal peak location estimator specifically for triangle waves in AGWN, compared with peak estimation for generic waveforms?

For a "smooth" waveform in noise, a common peak locator method might involve a linear-phase low-pass filter or an overdetermined polynomial curve fit (by regression, et.al.), followed by a maxima/minima search.  But if the signal is known to be a sharp triangle wave (which may imply that any sampling included/aliased some non-filtered frequency content above Fs/2, if that is relevant), or portion thereof, and not a low order polynomial or other smooth waveform, can this information be used to improve the triangle peak location estimation?

If so, how?

Added: Less than 1 period of the triangle wave may be available.  But assume at least one extrema is in the data window.

You may be able to determine this in the frequency domain. The noise signal is broad band and has low spectral density at any one frequency. The triangle, however, has all energy concentrated in just a few frequencies, so these should stick out in a short term Fourier spectrum. You can match against your expected spectrum (odd harmonics, amplitude proportional to 1/n^2).

The harmonics of a triangle are all in phase. Since your triangle wave is arbitrarily located inside your analysis window, the harmonics will appear to have a linear phase. You can probably model your measured phase (after proper unwrapping) as a linear phase using a simple least squares error approach. This should give you the exact position of the peak of the triangular wave with respect to the beginning for your analysis window.


SNR and noise in Simulations I've always taken the definition of noise and SNR for granted in doing BER simulations in an AWGN channel. I would usually use the AWGN noise function and make a loop with varying SNR values.
However, I am looking at other coding examples and I found this:

sig= 1/sqrt(2*SNR);
n1 = randn(NN,1);
noise1=sig*n1;
N0=2*(sig.^2);

Is this the definition embedded in the AWGN function in MATLAB?
If not, then what is it and why is sig defined this way and what is it in relation to N0?

The AWGN variance is:

sig= 1/sqrt(2*SNR);

Then, you define a random sequence:

n1 = randn(NN,1);

So the noise signal is:

noise1=sig*n1;

And the PSD is:

N0=2*(sig.^2);

Have a look to this example for more info.


Plotting DNA chromatogram trace data Sanger DNA sequencing produces a chromatogram trace which can be visualized with a number of programs, including FinchTV or ChromasLite. The raw data consists of co-ordinates for each of the four DNA bases (A, C, G and T). However, plots are displayed as smoothed peaks, as shown in the image previously. Many programs allow the x- and y-axis of the plot to be scaled larger or smaller to change the shape of the plot. What mathematical method is used to plot smooth curves such as these from a (small) number of raw data points?

Implementation

Assuming you already have a line drawing routine, you just need to supplement that with some kind of interpolation. The curves are created by drawing enough interpolated short lines to make the result look smooth. A good starting point would be to use an existing interpolation routine, like the ones given by Paul Bourke here.

I'll illustrate this using the cubic routines he supplies, since those are some of the simplest that will still give reasonable results. Here's the first one (translated to python) for reference:

def cubic(mu,y0,y1,y2,y3):
    mu2 = mu*mu
    a0 = y3 - y2 - y0 + y1
    a1 = y0 - y1 - a0
    a2 = y2 - y0
    a3 = y1
    return a0*mu*mu2 + a1*mu2 + a2*mu + a3

Each routine has a parameter mu which represents the fractional part of the index you want to interpolate. Depending on the routine, the other parameters will be some number of samples surrounding the index in question. In the cubic case, you'll need four samples. For example, if your data is y[n], and you want the value at 10.3, mu would be .3, and you'd pass in y[9], y[10], y[11], and y[12].

Then instead of drawing a single line with endpoints, say, $(10,y_{10}) \rightarrow (11,y_{11})$, you'd draw a bunch of shorter ones using the interpolated values (e.g. $(10,y_{10}) \rightarrow (10.1,\text{cubic}(.1,y_9,y_{10},y_{11},y_{12})) \dots $). Obviously those points would need to be scaled to the $x$ and $y$ dimensions of the image to be rendered.

Theory

Now, since the page/routine I referenced doesn't cite any sources, it's worth explaining where those cubic routines come from (and how they work). Both the one I reproduced above, and the very similar Catmull-Rom spline he mentions just below it, are two specific cases of using the following cubic convolution kernel:

$$
\psi(x) =
\begin{cases}
(\alpha + 2)|x|^3 - (\alpha + 3)|x|^2 + 1,  & \text{ if } 0 \le |x| \lt 1 \\
\alpha|x|^3 - 5\alpha|x|^2 + 8\alpha|x| - 4\alpha,  & \text{ if } 1 \le |x| \lt 2  \\
0, & \text{ if } 2 \le |x|
\end{cases}
$$

The routine listed above corresponds to a value of $\alpha = -1$, and the Catmull-Rom spline corresponds to $\alpha = -1/2$.  I won't go into too much detail about how the general form of the kernel is derived, but it involves various constraints like making sure that $\psi(x)$ is one at zero, and zero at all other integers.

This is what it looks like:

The two choices for the value of $\alpha$ come from attempts to match various aspects of the sinc function, the ideal reconstruction kernel. Setting $\alpha = -1$ makes the derivative of $\psi$ match the derivative of the sinc function at $x = 1$, and making it equal to $-1/2$ provides the best low-frequency approximation instead. By all accounts, a value of $\alpha = -1/2$ has much better properties overall, so it's probably the best value to use in practice. A much more extensive discussion can be found in the following paper, starting on page 328:

  Meijering, Erik.  “A Chronology of Interpolation: From Ancient
  Astronomy to Modern Signal and Image Processing.”  Proceedings of the
  IEEE. vol. 90, no. 3, pp. 319-42.  March 2002.

Insight

Now, just looking at that kernel in relation to the actual implementation of the interpolation code, it may not be clear how the two are related. Basically, the interpolation process can be thought of as adding together shifted copies of the kernel, which are scaled by the samples of the data, like so:

In fact, if you have an implementation of the kernel, you can use it directly to do the interpolation, as follows:

def kernel(x, a=-1.0):
    x = abs(x)
    if x >= 0.0 and x < 1.0:
        return (a + 2.0)*x**3.0 - (a + 3.0)*x**2.0 + 1
    elif x >= 1.0 and x < 2.0:
        return a*x**3.0 - 5.0*a*x**2.0 + 8.0*a*x - 4.0*a
    else:
        return 0.0

def cubic(mu,y0,y1,y2,y3):
    a = -1.0
    result  = y0 * kernel(mu + 1, a)
    result += y1 * kernel(mu, a)
    result += y2 * kernel(mu - 1, a)
    result += y3 * kernel(mu - 2, a)
    return result

However, it's much less computationally efficient to do it that way. As a bridge from the direct kernel approach to the more streamlined one above, consider that with a little bit of algebraic manipulation, the first implementation can be put in the following form:

def cubic(mu,y0,y1,y2,y3):
    mu2 = mu*mu
    mu3 = mu*mu2
    c0 = -mu3 + 2*mu2 - mu
    c1 =  mu3 - 2*mu2 + 1
    c2 = -mu3 + mu2 + mu
    c3 =  mu3 - mu2
    return c0*y0 + c1*y1 + c2*y2 + c3*y3

In this formulation, the c0...c3 values can be thought of as the coefficients of an FIR filter which is applied to the sample values. Now it's much easier to see how to derive the routine from the kernel.  Consider the kernel with $\alpha = -1$, like so:

$$
\psi(x) =
\begin{cases}
|x|^3 - 2|x|^2 + 1,  & \text{ if } 0 \le |x| \lt 1 \\
-|x|^3 + 5|x|^2 - 8|x| + 4,  & \text{ if } 1 \le |x| \lt 2 \\
0, & \text{ if } 2 \le |x|
\end{cases}
$$

Now evaluate that kernel symbolically at various shifted offsets, keeping in mind that mu ($\mu$) ranges from 0 to 1:

$$
\begin{eqnarray}
\psi(\mu + 1) &=& -(\mu + 1)^3 + 5(\mu + 1)^2 - 8(\mu + 1) + 4 &=& −\mu^3 + 2\mu^2 - \mu & & (c_0) \\
\psi(\mu) &=& \mu^3 - 2\mu^2 + 1 &=& \mu^3 - 2\mu^2 + 1 & & (c_1)\\
\psi(\mu - 1) &=& (1 - \mu)^3 - 2(1 - \mu)^2 + 1 &=& -\mu^3 + \mu^2 + \mu & & (c_2)\\
\psi(\mu - 2) &=& -(2 - \mu)^3 + 5(2 - \mu)^2 - 8(2 - \mu) + 4 &=& \mu^3 - \mu^2 & & (c_3)\\
\end{eqnarray}
$$

Note that $\mu - 1, \mu - 2$ get "flipped" to $1 - \mu, 2 - \mu$ respectively because of the absolute value on $x$ in the kernel definition.  Now we have the exact polynomials that are used in the "FIR version" of the interpolation routine. The evaluation of those polynomials can then be made more efficient through standard techniques (e.g. Horner's method). Similar things can be done with other kernels, and there are other ways of constructing efficient implementations as well (cf. Digital Audio Resampling Home Page).


Given a downscaled image, produce Gaussian blur of the original image Let's say I a perform convolution of an image with a Gaussian kernel of some size.

Is it possible to produce the exact same result from a downscaled version of the image? If so, what's its minimal size?

I'm not looking for an exact answer, some directions where I should start looking are enough. (I have some basic DSP and image processing knowledge but not much beyond that.)

Edit: A closely related question is - does Gaussian blur lose information, i.e. is it impossible to reconstruct the original from the blurred image?

I don't think you can do that.
You can probably get a similar result, in that you won't be able to tell by eye that they are different.

The reason I don't think you can do what you have in mind is this:
Imagine that you downscale (downsample, basically throw away pixels) an image, apply a blur, upscale the image.  Suppose now that you apply a blur to the original image.  If those two are the same, then you could deconvolve both images and from both of them recover the original image.

This would be the best compression format ever.  Downscale by 16X, blur.  Small file, send by email.  Upscale by 16X, deblur, original image.  Hell, why stop at 16X?  You could compress all of your images to a single pixel each and be able to recover them.

Unfortunately, things don't work that way.  Downscaling will (must) throw away data, and you can't recover it.

So, the downscaled and blurred image will be different from the blurred original.


Can I model process noise as a known "error" in my dynamics while designing a Kalman Filter? Consider I am modelling the dynamics of a robot and using a Kalman filter to obtain estimates of some state. I have certain terms in my equation which correspond to data not accessible to this robot ( states of other robots etc).

1) Is it fair to model these as process noise by assuming these terms to evolve based on some random process and act "Gaussian like" ? This also raises another important question:
2) Why does the Kalman Filter require us to have a positive definite covariance associated with the process noise. How do I interpret this in the real world when I am writing down process noise as some unmodeled physical terms ?

The Extended Kalman filter (EKF) effectively does that: Unmodelled nonlinearities are taken account of by assuming higher process or measurement noise, depending on where the nonlinearity appears.


Using spectral coherence for computing similarity between time series data I don't really know anything about signal processing. I need to find similairty between two 1-dimensional time series. I read somewhere that spectral coherence could be used to find similarity between two time series. But i am not able to get around how to use it. Can some one tell how i can use it to compute similarity between time series.

One way to check for similarity between time series is to use the cross-correlation:
$$
\rho_{xy}(\tau) = \frac{1}{\sigma_x \sigma_y} E[ (x[n] - \mu_x)(y[n] - \mu_y) ]
$$
where $\mu_{x,y}$ is the mean of $x$ or $y$ and $\sigma_{x,y}$ is the standard deviation of $x$ or $y$.

This gives us a value between -1 (completely anti-correlated) and +1 (completely correlated).

Sometimes, this doesn't work as well as we'd like because there is some filtering operation relating the two time series, so while there is a strong relationship between the two time series, the cross-correlation is not very close to 1.

As a result, the spectral coherence is sometimes used:
$$
C_{xy}(f) = \frac{|G_{xy}(f)|^2}{G_{xx}(f)G_{yy}(f)}
$$
where $G_{km}$ is the cross-spectral density between $k$ and $m$.

This is now a function of frequency, and the closer it is to 1 at a given frequency, the closer the two signals are related.

I've put together a simple example in R which illustrates the difference between cross-correlation and (spectral) coherence.

The top graph shows the two time series.  The middle graph shows the cross-correlation is less than 1.  The bottom graph shows the spectral coherence as a function of frequency.

You can see that the time series are not as similar in the higher frequencies as they are at low frequency.

Sometimes, the integral of these values is used over a limited bandwidth so that a single figure of similarity is arrived at.

R Code Below

  #15466

  x1 <- rnorm(1000)

  h <- c(1,3,2,1)/7
  x2 <- filter(x1,h, circular=TRUE) + rnorm(1000)/10

  par(mfrow=c(3,1))

  plot(x1,type="l",col="blue")
  lines(x2, col="green",lwd=5)
  title("Original and filtered noise signals")

  ccf(x1,x2)

  Ixy <- spectrum(cbind(x1,x2), plot=FALSE,spans=c(5,7))

  plot(Ixy$coh)h)


application of wavelet to vibration signal I ask you to help me to interpret this figure corresponding to a diagnosis of gear failure using CWT (Haar);
I have already worked on the analysis of the vibratory signals using the FFT; But since the signals from the bearings and the gears are non-stationary it is considered necessary to apply WT. So I started with simulated signals of gear failures to evaluate the sensitivity of this method compared to that of Fourrier

I've been working on this area for a while. The wavelet transform is a more powerful tool for time-frequency (TF) & feature extraction since the FFT or STFT always suffers from the resolution issues in TF domain. As a starter, I think the short videos provided in MATLAB is a very good point to start with. I paste the 1st part here and you can also find the rest parts from it.

https://www.mathworks.com/videos/understanding-wavelets-part-1-what-are-wavelets-121279.html

Bests,
GL


DFT-like transform using triangle waves instead of sin waves We know that DFT (discrete Fourier transform) breaks down a signal into multiple frequencies of sine waves. Does there exist a transform that does the same thing, but for triangle waves?

For my purposes, im only talking about 1-d signals (like voltages, etc). I'm studying historical stock market data, and I just want to look at reversals in certain stocks. In other words, I want to perform a "low-pass" on the stock price using this transform.

Edit:
If yes, how can I do it?

The closest orthogonal transform I know of that might meet your needs is the Slant Transform. It's based on sawtooth(ish) waves, but some of the basis functions do resemble triangle waves:

(source: Applied Fourier transform)

It was developed for image coding/compression, but it seems like a reasonable first approach for the analysis of long-term linear trends/reversals in financial data. It doesn't seem like many of the key papers describing the transform are available [for free] online, but the following paper probably has sufficient detail to implement something:

  A Truncation Method for Computing Slant Transforms with Applications
  to Image Processing. M. M. Anguh, R. R. Martin. IEEE Trans.
  Communications 43 (6), 2103-2110, 1995. (author link) (pdf link)

Specifically, see Section III which gives the recursion relations used to construct the transform matrix.


How to find the time instant where my sent signal and received signal can synchronize? I am having two .wav files, one is modulated wave that is made into a .wav file and another is a recieved .wav file which is recorded from microphone while playing first modulated .wav file. So now that I have two .wav files, one sent and one received. I am trying to synchronize the two files so that I can find the exact time instant where my received file gets in sync with sent file.

I am using, received_wave = thinkdsp.read_wave(filename = '.wav')
to read the wave from file and modulated_wave.write(filename = '.wav') to write the wave into file. I am using Think DSP module in python for doing these.

So this received_wave and modulated_wave has two attributes ys and ts where ys is amplitude and ts is time stamps. I want to find the ts such that both wave will synchronize.

for this right now iam using crosscorrealtion of modulated wave(with sync bits) and recieved wave(having both sync bits and data that iam sending such a way that start sync bits then data bits then again end sync bits) so when i crosscorrelate start sync modulated wave(modulated wave using start sync bits only)and recieved wave i can find the value in crosscorrelation matrix where it is almost equal to the auto correlation of sync_wave.

is this approach correct? or shud i check for maximum value of crosscorrelated matrix then find the index ?

This is usually addressed via looking for a known pattern in the data frame.

The modulation scheme is not mentioned but when synchronisation is required, a Phase Locked Loop (PLL) is most commonly employed. The PLL will track the received signal and tune a local oscillator to its frequency AND phase. When lock is achieved the signal that a local oscillator produces (at the side of the receiver) is perfectly "aligned" (or synchronised) with a component of the received signal (for example, the carrier). At this point, demodulation of actual information can start.

This still leaves us with the problem of knowing which bits in this stream of zeroes and ones produced by the demodulator are the actual message bits. For this reason, a known sequence of bits is used to mark the start of the actual frame. For more information, please see section "Preamble" from this link.

Of course, the PLL tracks the received signal continuously and counteracts phenomena such as Doppler (for example), which might have a significant effect in acoustic communications. So, in general, it would be worth implementing the PLL part.

Alternatively and if you are only interested in a play-back -> Record -> Demodulate kind of operation, then you can insert a sequence of symbols at the beginning of the transmitted information to force the modulator to produce a short segment of "known signal". A sequence of 10 zeroes for example, could be used to produce a segment of a clean sinusoid. At the receiving end and once you have recorded your signal, you can use cross-correlation to align (in an one-of manner) the known signal with the received and then use this (static) information to demodulate the rest of the received signal. This is still, in a very rudimentary way, what a PLL would do anyway.

Hope this helps.


difference between conventional analog butterworth filter and that implemented through IIR I am confused what is difference between conventional analog butterworth filter and that implemented through digital technique using IIR

If you have designed the IIR Butterworth filter using the bilinear transform method, then due to the logarithmic frequency mapping (warping), the frequency and phase responses will be distorted...


Digitally undoing a hardware filter I have a bunch of data that went through a hardware filter before being digitized. The data has been taken this way for quite some time, and everyone is fine with that. Now, there is a new application in development that requires unfiltered data. Rather than wait a year to collect new data, we would like to "unfilter" this old data. The transfer function of the hardware filter is: $$H(s)=\dfrac{24.2s+1}{2.2s+1}$$

My first throught was to take the inverse of this transfer function, and go straight to the z domain by applying the Eurler approximation. $$f_s=30Hz$$

$$H(z)=\dfrac{67-66z^{-1}}{727-726z^{-1}}$$

Converted to the difference equation:

$$V_{in}[n]=10.8507V_{out}[n]-10.8538V_{out}[n-1]+0.9851V_{in}[n-1]$$

I've run data through that difference equation in Scilab, and got a result. I want to validate the difference equation by running the result through the difference equation of the original filter, but I haven't been able to get the original data out of the second difference equation.

Which leads me to question, is what I've done here even valid? Is what I'm trying to do even possible?

Using the Euler approximation, the mapping between the $s$-plane and the $z$-plane would be:

$$
  s \rightarrow \frac{z - 1}{T}
$$
where $1/T = f_s = 30$ (thus $s \rightarrow 30(z-1)$).

Substituting back in the original transfer function yields:

$$
\begin{align}
  H(z) &= \frac{24.2\cdot 30 (z-1) + 1}{2.2\cdot 30 (z-1) + 1} \\
       &= \frac{726 (z-1) + 1}{66 (z-1) + 1} \\
       &= \frac{726z - 725}{66z-65} \\
       &= \frac{726 - 725z^{-1}}{66-65z^{-1}}
\end{align}
$$

Which converted to a difference equation is:

$$
V_{in}[n] \approx 0.0909V_{out}[n] - 0.0895V_{out}[n-1] + 0.9986V_{in}[n-1]
$$

As far as whether retrieving the original signal is possible partly depends on whether the original signal did not in fact include significant energy above $f_s/2$.


Confusion understanding causality? I already know the simple definition that causal system is the one that does not depend on future values of input but today i was confused as i came across a new definition of causality after reading "signal processing first "about causality  as shown in attached snapshots

I have attached two pages snapshots and i have drawn a thin red line to left side of common para, para that occurs in both snapshots

As shown in common para in  second line ,"cause does not precede corresponding effect"

This statement is about causal filter but what does this statement means in simple words?

Both snapshots contain a table ,upward snapshot contains table of non causal filter

Bottom snapshot contains table of causal filter

We can see from table of causal filter that both input and start at same value of time (n=0)while table of non_causal filter shows that output starts at earlier value of time(n=-2) while input starts at a later value of time(n=0),thus
Indicating that output(effect)precede input(cause) in case of non-causal filter

The statement is not correct. In a causal system, the cause (i.e. input) does indeed precede (come before in time) the corresponding effect (output).

The correct statement is the cause doesn't follow the corresponding effect.

"Cause creates the effect", so in a causal system, Cause ($x[n]$) will always come first and then the corresponding effect ($y[n]$).

  A causal $h[n]$ can not use samples $x[n+1], x[n+2], x[n+3], ...etc$ in order to create $y[n]$. $h[n]$ can only use $x[n], x[n-1], x[n-2],...etc$.


FFmpeg audio filter pipeline for speech enhancement Using ffmpeg, I'm recording a talk, its audio with a semi-distant microphone. It makes sense to somewhat enhance the speech clarity, reduce background static and environmental noise, etc., in the most generic terms.

What combination of ffmpeg audio filters am I best advised to use to enhance recorded speech?

I'm currently using

highpass-frequency=300,
lowpass-frequency=4000,

as that is the range speaking voice is most prominent in,

bass=frequency=100:gain=-50,
bandreject=frequency=200:width_type=h:width=200,

to additionally attenuate the bass frequencies, and

compand=attacks=.05:decays=.05:points=-90/-90 -70/-90 -15/-15 0/-10:soft-knee=6:volume=-70:gain=10

to apply dynamic range compression and supposedly make the softer passages (above -70dB) louder and the louder ones more even.

I understand digital filters are not magic, but I'm looking for something more from ffmpeg's (or LADSPA's) repertoire to enhance the sound more, make it more resilient to noise, sporadic coughs, perhaps lose the acquired metallic aftersound ...

The Linux command line I'm using to test:

filters='
...above filter lines, in sequence...
'

arecord -q -r 11025 -f s32_le -c 1 -t wav - |
    ffmpeg -f wav -i - -af "$filters" -f wav - |
    aplay

I just came across one article related to the compand command that might be useful (just in case someone looking for help sees this article).

From the article:

  To test this, and my understanding of compand, I added a very simple filter to remove the quiet section of the audio by decreasing their volume dramatically:

ffmpeg -i in.mp3  -filter_complex \
  "compand=attacks=0:points=-30/-900|-20/-20|0/0|20/20" \
  out.wav

  attacks=0 means that I wanted to measure absolute volume, not averaging the sound over a short (or long period of time). When the speaker suddenly yells, or talks back and forth with someone off mic, I want the volume adjustment to be immediate. The downside is that you may hear the volume being clamped.

  points is the actual volume mapping function, and I’ll walk through it:

  I added a mapping of -30/-900, which means that volume below -30db in the original input track gets converted to -900db (completely silent).


GCC-PHAT (Generalized cross correlation MATLAB) I've been trying to implement generalized cross correlation with a PHAT weighting function for a while now, and cannot get it to work.  I have tried correlating using MATLAB's xcorr.m file, and it works for getting a correct delay tau (on simulated sinusoidal signals).

The code for the signals are thus:

Fs = 8000;
dt = 1/Fs;%0.125e-3
f1 = 100;
tdelay = 0.625e-03;%try different values
t3 = (0:dt:(1)-dt)';
x3 = cos(2*pi*f1*t3);
x4 = cos(2*pi*f1*(t3-tdelay));

As can be seen in part of the source code for xcorr.m, the cross correlation is implemented thus:

%Transform both vectors

X = fft(x,2^nextpow2(2*M-1));
Y = fft(y,2^nextpow2(2*M-1));

% Compute cross-correlation

c = ifft(X.*conj(Y));

According to definitions of GCC-PHAT, the only addition I needed to make was to divide the product by its own magnitude, before taking the ifft.  Here is my version with this change.

%Transform both vectors
X = fft(x,2^nextpow2(2*M-1));
Y = fft(y,2^nextpow2(2*M-1));

% Compute cross-correlation

R = X.*conj(Y);
c = ifft(R./abs(R));

However I always end up with a tau of zero with the PHAT weighting!  On closer examination of the array produced as a result of this division, it seems as if the first value of R is a real value (with no imaginary component), and so when divided by its magnitude it becomes 1.  All the other values in array R are complex, so do not end up as 1 when divided by their own magnitude and so end up with a value of <1.

This can be seen below, for the 1st 10 values of R.

    K>> R(1:10,1)

ans =

  0.000000000000000 + 0.000000000000000i
 -0.494299608718696 - 0.003002230689022i
 -0.002678647083223 - 0.000032538742345i
 -0.488954228290329 - 0.008909374553649i
 -0.010656518992354 - 0.000258902698589i
 -0.478379290671260 - 0.014528074329782i
 -0.023760667475633 - 0.000865926459320i
 -0.462803929640386 - 0.019677623220519i
 -0.041707017319469 - 0.002026674993917i
 -0.442565618721743 - 0.024194329448597i

K>> abs(R(1:10,1))

ans =

   0.000000000000000
   0.494308725968464
   0.002678844707371
   0.489035391682370
   0.010659663580139
   0.478599844010494
   0.023776441018801
   0.463222070011989
   0.041756229537848
   0.443226457301486

K>> R(1:10,1)./abs(R(1:10,1))

ans =

  1.000000000000000 + 0.000000000000000i
 -0.999981555555690 - 0.006073594357736i
 -0.999926227844417 - 0.012146557900713i
 -0.999834033705084 - 0.018218261306199i
 -0.999705001216859 - 0.024288074069393i
 -0.999539169638280 - 0.030355367874845i
 -0.999336589392997 - 0.036419515378054i
 -0.999097322000240 - 0.042479891383435i
 -0.998821440083941 - 0.048535871565711i
 -0.998509027227826 - 0.054586834901284i

As can be seen above, the largest value ends up being at the 1st index, when we divide R by its own magnitude.  So once the ifft is taken, the highest value is ALWAYS at the beginning of the array, which gives a lag and timediff of zero...even when I set the delay between the 2 identical signals at the beginning to a nonzero value (e.g. delay = 0.75e-03).

What am I doing wrong?? Any help appreciated.

Rory

P.S. if anyone is wondering why I'm bothering with the PHAT weighting, its because it should give far better results (in theory) in a real life scenario, for TDOA.

I remember I was facing a similar problem some time back. I am pretty sure it has to do with the DC and nyquist values when performing the equalisation and then the inverse. The first element in the array which is becoming 1.0 is the DC value.

I have this code which I’m pretty sure works for GCC-PHAT:

function [r,tau] = gccphat( x, y, fs )
M = max(numel(x),numel(y));

%%Transform both vectors
% X = fft(x,2^nextpow2(2*M-1));
% Y = fft(y,2^nextpow2(2*M-1));
%
% % Compute cross-correlation
%
% R = X.*conj(Y);
% c = ifft(R./abs(R));

%%
N = 2*M-1;
Nfft = 2^nextpow2(N);

R = bsxfun(@times, ...
        fft(y,Nfft), ...
        conj(fft(x,Nfft)));
rtmp = fftshift( ...
        ifft(exp(1i*angle(R))) ,1);
r = rtmp(Nfft/2+1-(M-1)/2:Nfft/2+1+(M-1)/2,:);

lags = (-(N-1)/2:(N-1)/2).';
lags = lags/fs;
[~,idx] = max(abs(r));
tau = N/(2*fs)+lags(idx);

end


How did Scientists Capture First Images of Molecules Before and After Reaction? "Using a state-of-the-art atomic force microscope, the scientists have taken the first atom-by-atom pictures, including images of the chemical bonds between atoms, clearly depicting how a molecule's structure changed during a reaction. Until now, scientists have only been able to infer this type of information from spectroscopic analysis."

http://www.sciencedaily.com/releases/2013/05/130530142007.htm

Are there any papers how this thing actually works?

The seminal paper you're looking for is:

Binnig, G., & Quate, C. F. (1986). Atomic force microscope. Physical Review Letters, 56(9),    930-933.

Here's the short version, along with some speculation since I don't have their manuscript.

Atomic force microscopy (AFM) is a scanning probe microscopy technique, meaning a sharp probe is scanned in a raster pattern over/in contact with the sample. A probe my look like this.

I will refer to the point as the tip, and that is where the carbon monoxide molecule mentioned in the article resides. The substrate the tip is attached to will be referred to as the cantilever.

A laser is bounced off the top of the cantilever, onto a dicrotic mirror, and finally into a photodetector, as shown below.

The researcher has control of the laser and mirror alignment. Alignment procedure is to place the laser centered, and as close to the end of the cantilever as possible, for the maximum deflection signal. The mirror is then adjusted to put the reflected laser on the most sensitive part of the photodetector.

The article mentions using non-contact AFM. That is about as meaningful as saying "we can the signal through a digital filter." There are multiple non-contact modes, but my guess (without having access to the manuscript) is they used a form of tapping mode.

The cantilever has a self resonant frequency, typically around 300 kHz. The cantilever is driven to resonance, and scanned over the sample. The deflection signal generated is an amplitude modulated version of the topography. It is also possible to obtain information about material properties this way.

That is the short version. If I get access to the manuscript or someone asks questions, I would be happy to add more information about how this system works.


Why am I not getting the intended coefficients in this 2D Fourier demonstration? I am trying to demonstrate how the 2D Fourier decomposition of an image works with Matlab and a very, very simple example.

I create a 4 x 4 pixel image with cosine basis vectors as such:

%Basis vectors with coefficients:

f = @(x, y,a0,a1,a2,a3,b0,b1,b2,b3,c0,c1,c2,c3,d0,d1,d2,d3)...
a0*cos(0*x + 0*y) + a1*cos(0*x + 1*y) +...
a2*cos(0*x + 2*y) + a3*cos(0*x + 3*y) + ...
b0*cos(1*x + 0*y) + b1*cos(1*x + 1*y) + ...
b2*cos(1*x + 2*y) + b3*cos(1*x + 3*y) + ...
c0*cos(2*x + 1*y) + c1*cos(2*x + 1*y) + ...
c2*cos(2*x + 2*y) + c3*cos(2*x + 3*y) + ...
d0*cos(3*x + 1*y) + d1*cos(3*x + 1*y) + ...
d2*cos(3*x + 2*y) + d3*cos(3*x + 3*y)

%Matrix of image space:
[X,Y] = meshgrid(0:3);

%Creating the image from Fourier space by assigning coefficients:

F = f(X,Y,1,1,-1,1,1,1,1,6,3,1,1,1,7,1,1,1);
Fourier = fft2(F);

F =
 27 -5.53409 3.93452 -10.2602
 -0.177994 -12.1365 8.43196 1.46219
 3.50819 6.86722 2.91194 -3.6156
 -18.1397 2.3525 -6.23362 5.35198

Fourier =

 (5.72275,0) (3.14567,1.38924) (36.7478,0) (3.14567,-1.38924)
 (5.46847,-14.2485) (33.0684,2.46052) (43.5603,-51.0061) (11.87,-9.05285)
 (43.9012,0) (44.1778,-31.8071) (63.0469,0) (44.1778,31.8071)
 (5.46847,14.2485) (11.87,9.05285) (43.5603,51.0061) (33.0684,-2.46052)

The image is:

Now I naively want to get back the coefficients with Fourier = fft2(F), but it is not happening - to begin with, there are imaginary components, which I presume are part of a phase shift I did not introduce (I only used cosines in the creation of the original image).

What are my mistakes and misconceptions undermining this demonstration?

I see that the basis vectors can actually be expressed as complex numbers, e.g.

$$\cos(3x + 1x) = \frac 1 2 e^{-3 i x - i y} + \frac 1 2 e^{3 i x + i y}.$$

Yet, predictably, expressing the function with complex terms does not change much - below is the function expressed with complex terms for access:

f = @(x, y,a1,a2,a3,a4,b1,b2,b3,b4,c1,c2,c3,c4,d1,d2,d3,d4)...
a1*(1/2*exp(-0*i*x-0*i*y) + 1/2*exp(0*i*x+0*i*y)) + a2*(1/2*exp(-0*i*x-1*i*y) + 1/2*exp(0*i*x+1*i*y)) +...
a3*(1/2*exp(-0*i*x-2*i*y) + 1/2*exp(0*i*x+2*i*y)) + a4*(1/2*exp(-0*i*x-3*i*y) + 1/2*exp(0*i*x+3*i*y)) + ...
b1*(1/2*exp(-1*i*x-0*i*y) + 1/2*exp(1*i*x+0*i*y)) + b2*(1/2*exp(-1*i*x-1*i*y) + 1/2*exp(1*i*x+1*i*y)) + ...
b3*(1/2*exp(-1*i*x-2*i*y) + 1/2*exp(1*i*x+2*i*y)) + b4*(1/2*exp(-1*i*x-3*i*y) + 1/2*exp(1*i*x+3*i*y)) + ...
c1*(1/2*exp(-2*i*x-0*i*y) + 1/2*exp(2*i*x+0*i*y)) + c2*(1/2*exp(-2*i*x-1*i*y) + 1/2*exp(2*i*x+1*i*y)) + ...
c3*(1/2*exp(-2*i*x-2*i*y) + 1/2*exp(2*i*x+2*i*y)) + c4*(1/2*exp(-2*i*x-3*i*y) + 1/2*exp(2*i*x+3*i*y)) + ...
d1*(1/2*exp(-3*i*x-0*i*y) + 1/2*exp(3*i*x+0*i*y)) + d2*(1/2*exp(-3*i*x-1*i*y) + 1/2*exp(3*i*x+1*i*y)) + ...
d3*(1/2*exp(-3*i*x-2*i*y) + 1/2*exp(3*i*x+2*i*y)) + d4*(1/2*exp(-3*i*x-3*i*y) + 1/2*exp(3*i*x+3*i*y))

If I run imshow(ifft2(Fourier)) I do get the original image with the initiall 4 x 4 setup.

So I assume I am on the right track, but I can't match the coefficients in the 2D DFT.

I have been able to do it for the DC component:

g = @(x, y, a1)...
dot(F(:), (a0* cos(0*x + 0*y))(:))
g(X,Y,1)
% ans = 5.72275

Which performs the Frobenius product of the image matrix multiplied by the matrix of the DC component:

$$a_0 \cos(x + y).$$

The result is 5.72275.

But I haven't been able to do the same with, say, the coefficient for

$$b_3\cos(x +3y)$$

with

q = @(x, y, b3)...
dot(F(:), (b3* cos(1*x + 3*y))(:))
q(X,Y,6)
% ans = 449.363

which yields 449.363, while there are no coefficients greater than 80 in the Fourier matrix Fourier. Evidently I am messing up treating the complex values correctly...

Your misconception is to assume that you can compute Fourier coefficients with the fft. You cannot (really). Fourier coefficients are found by integrating a continuous time (periodic) function after multiplying it with a harmonic kernel. The cosine function gives two real-valued Fourier coefficients equal to half its amplitude in this case.

You are considering a matrix though which means that you have sampled your cosine function at a discrete grid (which is not a problem in your case) and truncated it to a finite width (which is exactly the problem you have). Then you are computing the FFT which is a fast way of computing a DFT (Discrete Fourier Transform). The DFT (of discrete vectors) is not the same as a Fourier transform (which is defined on continuous functions).

That said, your DFT can approximate Fourier coefficients, as the Fourier integral can be approximated by a sum, which is precisely the definition of the DFT then. The approximation errors that happen in this case are sampling and truncation.

For functions that are strictly bandlimited (like your cosine), sampling does not cause any error as long as you satisfy the Nyquist criterion. Conversely, for functions that are strictly time-limited, truncation does not cause any error as long as the whole function fits in your window. Unfortunately, it can be shown that no function can be time-limited and band-limited at the same time. For this reason, the DFT can never accurately predict the result of the Fourier transform.

There is hope though, since you're not after a Fourier transform but instead you're after Fourier coefficients. These are defined on periodic continuous-time functions. For these, the DFT can give you accurate predictions of the Fourier coefficients, as long as you make sure that the samples you are considering contain a whole period of the signal. This needs to be accurate, i.e., if you periodify your extracted window it should look like the original function. For instance, if the period of your function is 1 and you want to extract four samples, they need to be at 0, 0.25, 0.5, 0.75 (i.e., not 0, 1/3, 2/3, 1, which is a very common mistake in this setting).

Your Fourier coefficients will then be real-valued if and only if your function is symmetric, which in the sampled case means x[n] = x[N-n] for n=1, 2, ..., N-1.  If you sample your cosines correctly, this should work.


Anisotropic Diffusion Filter - Intuition Behind Parameters I need to know what is iteration and divergence in anisotropic diffusion filter technique.

  Isotropic diffusion $$\frac{\partial I(x, y, z)}{\partial t}={\rm div}\left[c\cdot \nabla I\left(x, y, z\right)\right], \quad
\text{where } c \text{ is the diffusion coefficient}$$
  Anisotropic diffusion $$\frac{\partial I(x, y, z)}{\partial t}={\rm div}\left[g\left(\left\| I\left(x, y, z\right)\right\|\right)\cdot \nabla I\left(x, y, z\right)\right], \quad \text{where } g \text{ is the anisotropic diffusion coefficient }\textbf{(Edge stopping function)}$$

Here $t$ refers to iteration, what is this iteration ? How is it related to filtering ? I know that iteration refers to number of rounds but how it is related with filtering ?
Also explain to me why we use divergence in these equation what is the purpose of divergence.

This paper - The Structure of Images describes quite nicely to get to the diffusion equations from the point of view of a filter. It's about the isotropic case, but anisotropic is really just an extension of that. For anisotropic diffusion see Scale Space and Edge Detection Using Anisotropic Diffusion as Deniz suggests this one is quite good as it also includes discretised versions of the equations.

The basic idea is to solve the heat equation where $t$ is the scale of the image, i.e. the amount of filtering to apply. For anisotropic diffusion the diffusion coefficient $c$ is replaced by $g$ which is dependent on the image gradient and therefore filters preferentially depending on the gradient.

The divergence here really just arises from the definition of the heat equation.


What is the low-pass filter that gives minimal transient response? As title says, what is the low-pass filter that gives minimal transient response? By minimal transient response, I mean minimal settling time with reference to very tiny amount of divergence of steady-state values.

I am asking for analog case, but digital case would also be fine.

As mentioned in the comments there is no unique answer to your question, because you're asking for a contradiction. A good low pass filter (i.e. one with a fast transition between pass band and stop band, and with high stopband attenuation) will exhibit a considerable transient response. On the other hand, a filter with a very smooth step response is worse as a frequency selective filter than other filters of the same order. However, a Bessel filter might be a good compromise between these two conflicting goals. It performs reasonably well as a low pass filter, and it has very little overshoot in its step response.


Difference between double-sided and single-sided AWGN noise after bandpass filtering? Consider a bandpass signal $s(t)$ with bandwidth $W$.
After bandpass filtering, let the output signal be $r(t)=s(t)+n(t)$
I have read a paper that denotes $n(t)$ as Gaussian noise with one-sided power spectral density $N_0$. Therefore, the noise power is $\sigma^2=\mathbb{E}\{n^2(t)\}=N_0W$.
What would be the purpose of denoting the noise as single-sided?
It seems that if we consider the noise as double-sided with power spectral density $\frac{N_0}{2}$, noise power is still $\sigma^2=\mathbb{E}\{n^2(t)\}=N_0W$ since we have to integrate over the negative frequencies and the positive frequencies.
What is the purpose of describing AWGN noise as single-sided versus double-sided? Considering the case of real signals, do both end up giving the same results?

Both descriptions give the same result. We often use the one-sided noise power spectral density (PSD) because for real-valued processes the negative frequencies are redundant, so defining the PSD for positive frequencies is sufficient. You just have to scale the noise power spectrum such that integrating the one-sided power spectral density (PSD) over the positive frequencies gives the same result as integrating the two-sided PSD over positive and negative frequencies. I.e., for white noise with a constant PSD, defining the one-sided PSD as $N_0$ means that we would need to define the two-sided PSD as $N_0/2$ such that the noise power remains the same: $N_0W=N_0/2\cdot 2W$


Determine the optimum receiver and the corresponding $P_{eM}$ for an AWGN channel I have a source that emits $M$ equiprobable messages, which are assigned signals $s_1, \dots,s_M,$ that are equidistant by $a$. That is, if we plot the $s_k$ signals in a horizontal axis they are dots separate by the same distance $a$. And so I would like to determine the optimum receiver and the corresponding error probability $P_{eM}$ for and AWGN channel as a function of $E_b/\cal N$.

As has been referenced in the comments already, you're describing a pulse-amplitude modulation (PAM) signal constellation. The problem, as you've framed it, seems to suggest the AWGN vector channel (where symbols are described using discrete values $s_1$, $s_2$, and so on), in contrast to the waveform channel, where symbols are expressed using waveforms that are transmitted over some period of continuous time.

While in the real world, the waveform channel is often used as a communication medium, analysis of performance of digital communication systems is typically related back to the vector channel for ease of analysis. This document has a discussion of the distinction between the two channels and how the conversion between the two is modeled mathematically.

With that said, you can frame the problem rather simply. Assume that at each symbol time, your receiver observes a value $x$, which is defined as follows:

$$
x = s_i + n,
$$

where the value of $i$ indicates the symbol that was transmitted. $s_i$ refers to the value that would be expected at the receiver in the case of no noise (i.e. $s_i$ refers to the amplitude that is transmitted for symbol $i$). I will state without proof that the optimal decision rule (for minimum probability of bit error) for the case of equiprobable input symbols and white Gaussian noise is to pick the symbol value whose constellation position is closest to the observed point. This result follows from Bayesian decision theory.

The exact expression for the bit error rate will depend upon the coding of bits onto the symbols in your constellation. Borrowing a figure from an above comment:

These are examples of Gray coded constellations, such that adjacent symbol positions only differ by one bit. These codings are used to minimize the probability of bit error for multi-bit-per-symbol modulations. You might ask: why do you need to know the symbol encoding in order to calculate the system's error performance?

The answer might not be obvious: the overall bit error probability is a single number that is averaged across the ensemble of all possible transmitted symbols and observed channel noise. However, it should be clear from the above figure that the bit error probability must be different based upon which symbol is actually transmitted.

  Example: Assume the $M=8$ case illustrated above, with a transmitted symbol of $000$.

  The noise observed at the receiver $n$ can take on any value $< 0$ without having any chance of causing a bit error! $n$ could be $-10000000a$ (which is likely not very probable), but $000$ would still be the closest constellation position and therefore the receiver would make the right decision. We clearly couldn't say this for another transmitted symbol like $110$, however, so the bit error probability must be dependent upon the symbol that was actually transmitted.

Assuming a particular coding of bits onto the constellation (e.g. the illustrated Gray code), you can construct the error probability conditioned on the actually-transmitted symbol. I will tackle the $M=4$ case from above because it is simpler than the 8-point constellation, but the same principle applies. While the spacing between points isn't shown on the graph (because I'm too lazy to draw one), the points from left to right should be located at $-1.5a$, $-0.5a$, $0.5a$, and $1.5a$, respectively. Remember that the optimum receiver chooses the constellation point that is closest to the observed value.

$$
P(\text{error | 00 transmitted}) = \frac{1}{2} \\
(P(\text{choose 01 | 00 transmitted}) + \\
2 P(\text{choose 11 | 00 transmitted}) + \\
P(\text{choose 10 | 00 transmitted}) )
$$

So we've very straightforwardly stated what the probability of a bit error is for this case. If $00$ was transmitted and $01$ is chosen by the receiver, then you have one bit error; if $11$ is chosen by the receiver, you have 2 bit errors; and so on. Since we know the decision rule employed by the receiver, we can state this more directly in terms of the observation. Remember that if the symbol $00$ was transmitted, then the observation should have the form:

$$
x = s_i + n = -1.5a + n
$$

$$
P(\text{error | 00 transmitted}) = \frac{1}{2} \\
(P(-a < x < 0\ |\ x = -1.5a + n) + \\
2P(0 < x < a\ |\ x = -1.5a + n) + \\
P(x > a\ |\ x = -1.5a + n) )
$$

We can rearrange the above to be solely in terms of the noise term $n$:

$$
P(\text{error | 00 transmitted}) = \frac{1}{2} \\
\left(P(0.5a < n < 1.5a) + \\
2P(1.5a < n < 2.5a) + \\
P(x > 2.5a)\right)
$$

On the AWGN channel, $n$ is a sample from a zero-mean white Gaussian random process. Assuming the noise process has variance $\sigma^2$, the above probabilities can be written in terms of the Gaussian distribution's tail probability function, known as the Q function:

$$
P(\text{error | 00 transmitted}) = \frac{1}{2} \left(\left(Q\left(\frac{0.5a}{\sigma}\right) - Q\left(\frac{1.5a}{\sigma}\right)\right) + 2\left(Q\left(\frac{1.5a}{\sigma}\right) - Q\left(\frac{2.5a}{\sigma}\right)\right) + Q\left(\frac{2.5a}{\sigma}\right)\right)
$$

In most cases, the error probability will be dominated by the closest constellation points (the noise is unlikely to be large enough to cause you to skip multiple points over in the constellation, so the corresponding probabilities will be small compared to those for the nearby points), but for strict correctness, you would need to include all of the terms. You can then repeat this analysis for each constellation point, calculating what the probability of making all of the incorrect choices are. Then, to arrive at a single bit error probability number, you use the law of total probability:

$$
P_e = \sum_{i=1}^{M} P(\text{error | }s_i\text{ transmitted}) P(s_i\text{ transmitted})
$$

Since all of the symbols are equally probable,

$$
P(s_i\text{ transmitted}) = \frac{1}{M}
$$

$$
P_e = \frac{1}{M} \sum_{i=1}^{M} P(\text{error | }s_i\text{ transmitted})
$$

And that is how you would calculate the bit error performance for an arbitrary PAM constellation. Note that it is inherently a function of the signal-to-noise ratio, implied via the proxy values of $a$ (which is related to the received amplitudes, which is related to the received "energy per bit" $E_b$) and the noise variance $\sigma^2$ (which is also often dubbed $N_0$ in digital communications analysis). You could manipulate the above for given values of $M$, $a$, and $\sigma^2$ to arrive at a final result that is a function of the typical metric $\frac{E_b}{N_0}$. That is left as an exercise to the reader.