Filtering a signal using Autoregressive (AR) filter and finding the coeff of AR filter using Yule Walker equation in MATLAB

Before I go into the specific issues, I'd like to note that there is a minor confusion with respect to M not representing the AR process order. More specifically an AR process of order M=4 would use the M previous values of y1 rather than M+1 values (in the loop with index j). So, with M=4 your implementation is really that of an AR process of order $p = 5$.

Coefficients and stability

  The problem is that I need to multiply the previous output y1 and current input x with a=0.1, otherwise it gives an unstable filter which leads to infinity. Please clarify which mistake I'm making while filtering x using AR filter or I'm doing it the wrong way?

The stability of the AR process can be determined from the location of the poles of the transfer function. To find those poles, we start off with the recurrence equation of the AR(p) process:

$$
  y[n] = a x[n] + a \sum_{j=1}^{p} y[n-j]
$$

From which the transfer function in the $z$-plane is:

$$
\begin{align}
  H(z) &= \frac{a}{1-a\sum_{j=1}^{p} z^{-j}} \\
       &= \frac{az^{p}}{z^{p}-a\sum_{j=0}^{p-1} z^j}
\end{align}
$$
Hence the poles are given as the roots of the $z^{p}-a\sum_{j=0}^{p-1} z^j$ polynomial. In your specific case of $p = 5$, the pole location are given by the following plot:

As you may notice, there is a pole near $z=2$ (more precisely at approximately 1.96595) when $a=1$ which indicates that the system is not stable. On the other hand, when $a=0.1$, all the poles are located inside the unit circle and the system is stable. Beware that even with $a=0.1$, starting with $p=11$, the system also has a pole outside the unit-circle and starts to be unstable.

In other words, there are no specific mistaking in your filtering code (ie. it correctly implements an unstable filter when a=1).

Yule-Walker equation

  But its not giving correct coefficients (or filter taps). Please somebody help me if there's something I'm missing in understanding the question fully or making some mistake somewhere in the code?

There are indeed a few issues with that part of the code. Perhaps the most significant issue is that the solution

$$
\begin{align}
    \vec{a} &= \mathbf R^{-1} \vec r \\
      &= \begin{bmatrix}
           r_0 & r_1 &\cdots & r_{p-1} \\
           r_1 & r_0 &\cdots & r_{p-2} \\
           \vdots & \vdots & \ddots & \vdots \\
           r_{p-1} & r_{p-2} &\cdots & r_0 \\
         \end{bmatrix}^{-1}
         \begin{bmatrix}
           r_1 \\
           r_2 \\
           \vdots \\
           r_p
         \end{bmatrix} \tag{1}
\end{align}
$$
where
$$
\begin{align}
  r_j &= \frac{c_j}{c_0} \\
      &= \frac{E\left\{y_n y_{n-j}\right\}}{E\left\{y^2_n\right\}}
\end{align}
$$

applies to an AR process of the form
$$
  y[n] = x[n] + \sum_{j=1}^p a_j y[n-j]
$$
where $x[n]$ is a zero-mean random variable. Unfortunately, the line

x = rand(1,1000);

Generates a uniformly distributed random variable in the range [0,1] with a mean of $\frac{1}{2}$. An example of random variable with zero-mean could be a uniformly distributed random variable in the range [-1,1]:

x = 2*rand(1,1000)-1;

For an input which does not have a mean of zero, or equivalently for an AR process of the form:
$$
  y[n] = K + x'[n] + \sum_{j=1}^p a_j y[n-j]
$$
where $x'[n] = x[n] - K$ is zero-mean, the equations can be updated to:
$$
\begin{align}
  \begin{bmatrix}
\mu \\
c_1 \\
c_2 \\
\vdots \\
c_p
  \end{bmatrix} &=
\begin{bmatrix}
  1   & \mu & \mu &\cdots & \mu \\
  \mu & c_0 & c_1 &\cdots & c_{p-1} \\
  \mu & c_1 & c_0 &\cdots & c_{p-2} \\
  \vdots & \vdots & \vdots & \ddots & \vdots \\
  \mu & c_{p-1} & c_{p-2} &\cdots & c_0 \\
\end{bmatrix}
\begin{bmatrix}
  K \\
  a_1 \\
  a_2 \\
  \vdots \\
  a_p
\end{bmatrix} \tag{2}
\end{align}
$$
where $\mu = E\left\{y_n\right\}$.

Implementation for zero-mean AR process

Finally there are a few small issues in the implementation of the equation. For one thing, the vector returned by autocorr(y2,p) contains p+1 values: from 0 to p lags. You need the first p values for the $\mathbf R$ matrix (leaving out lag p), and the last p values for the $\vec r$ vector (leaving out the zero lag value since $\vec r$ starts at $r_1$, not $r_0$).

So you can construct the Yule Walker equation with:

R = toeplitz(r(1:end-1)); % r(0), ..., r(p-1)
w = inv(R)*r(2:end);      % as Yule Walker eq: w=R^-1*r

So putting it all together, you should get something like this:

x = 2*rand(1,1000)-1; % uniformly distributed over [-1,1]
a=0.1
p=5;
M=p-1;
for i = 1:length(x)
  sum=0;
  y1(i)=0;
  for j = 1:M+1
    if (i-j>0)
      sum = sum + a*y1(i-j);
    end
  end
  y1(i)=sum +a*x(i);
end

y2 = awgn(y1,25);

r = autocorr(y2,p)';
% Or if you don't have autocorr, replace with the following 2 lines:
% r = xcorr(y2,p)';
% r = r(p+1:end);           % keep only non-negative lags

R = toeplitz(r(1:end-1)); % r(0), ..., r(p-1)

w = inv(R)*r(2:end);      % as Yule Walker eq: w=R^-1*r

Implementation for non-zero-mean AR process

The solution of equation (2) can be constructed in a similar fashion as for equation (1). The main thing to note is that in this case the $\mathbf R$ matrix elements $c_0, \dots, c_{p-1}$ are not normalized by $c_0$ (which would have yielded $r_0, \dots, r_{p-1}$). In this case you may replace autocorr with xcorr, which returns a sum of product without normalization (you would then need some adjustment since is xcorr is double sided, meaning it contains 2*p-1 values, including the values for positive, zero and negative lags). This would give you the following:

x = rand(1,1000);
...

mu = mean(y2);
r = xcorr(y2,p)';
r = r(p+1:end);           % keep only non-negative lags
R = toeplitz(r(1:end-1)); % r(0), ..., r(p-1)

A  = [1 mu*ones(1,p);mu*ones(p,1) R];
b  = [mu;r(2:end)];
w  = inv(A)*b;

% separate the bias term K from the AR process weights:
K  = w(1);
w  = w(2:end);


Algorithm to convert a GeoTIFF file to PNG image format

I can't put tiffs directly into a TFRecords format...

You can. TFRecords does not know what an "image" is, it simply stores a string of bytes. What this means to higher level applications is up to them.

When I imported a sample GeoTIFF file into python using rasterio I found that the range for pixel values was between [282, 1560], which seemed rather high.

Depending on the kind of data they contain, GeoTIFF images can have higher colour depth per pixel. For example, the SRTM data provide a digital elevation model of the Earth. Each pixel value is associated with the height of a "square" that is 30 by 30 meters (along the equator). The variation of that height over the whole planet is much larger than 256. SRTM's GeoTIFFs have a 16bit colour depth.

Now PNG or JPEG files seem to have pixel ranges of between [0, 255] or [0,1], from what I gather.

PNG and JPEG can save 16bit per pixel colour depth images and you can still work with them with libraries such as osgeo or GIMP.
Depending on the content of the GeoTIFFs though, the data in PNG format (lossless encoding) might not compress as much as you might expect and with JPEG (lossy encoding) you are going to have to tune the compression to avoid losing any necessary features. In the end, it might be better to just leave the data in GeoTIFF or read them in a numpy array and then serialise that straight into TFRecords.

This seems like a simple rescaling of the values to my naive eyes.
So my question is whether anyone could provide the basic algorithm to convert TIFF file pixels into PNG or JPEG pixels. I don't need to worry about the spatial coordinates in this case, so I can treat this as a straight TIFF --> PNG or JPEG conversion.

Indeed, it has to be a simple linear scaling but this is not a good idea for a number of reasons.
If the absolute values are of importance (for example) then the 8bit dynamic range is simply not adequate.
If your subsequent processing is about learning gradients (for another example), then you are going to have to normalise to the global extrema (not the maximum / minimum of each individual file) so that the relative slopes do not change. The danger in doing so is again dynamic range. Although it is possible to squash any numerical range in the interval $0 \dots 1$, when this interval is quantised in 8 bits it is only divided into 256 "steps". If two pixel values differ by less than $\frac{1}{256}$, then they will be assigned to the same "bin". Also, if the dataset contains "spikes", that is the majority of the values are between 20-100 (for example) but there is one value at 10000, then these spikes will "kill" the contrast when you normalise at 8bits. What you will get is a mostly black image with 1 white pixel. Again, it might be possible to remove those spikes but if you have to work with them (e.g. for a multispectral application), then normalising to an 8bit colour depth is not a good idea.
The other thing you would be losing if you went with this option is the transform data. That is, the transform that converts coordinates from image (pixel positions and values) to world coordinates (longitude, latitude, height).
Suppose that your GeoTIFF image ($I$) describes a digital elevation model that has both negative and positive values. The total range of the GeoTIFF would be $max(I)-min(I)$. If $min(I)$ is not zero then your $0$ point in the converted data will be "floating". It might be 20 on one image, 42 on another 100 on yet another and so on.
The transformation between two ranges is straigthforward: $O_{x,y} = O_{min} + \frac{I_{x,y}-min(I)}{max(I)+min(I)} \times range(O)$, where $O$ is the output image, $O_{min}$ is the new desired minimum and $range(O)$ is the desired range.
While the transformation is straightforward, the impact it has on the data and subsequent models can be severe, therefore my recommendation would be to preserve the original numbers as much as possible and only intervene if you are absolutely sure that the transformation will not affect your objectives.
Hope this helps.


How to upscale an image without affecting pixel aspect ratio?

Popular cheap methods are:

Letterboxing

Cropping

A combined method: cropping, overlapped by the letterboxed image

The combined method works surprisingly well in videos, and people are quite accustomed to viewing such videos.

Photoshop and others also have smart content-aware scaling/resizing algorithms that re-organize or generate image content to the new image size. I don't have PhotoShop to demonstrate it with, but here is a video from Adobe.


Indices in 2D Convolution Operations - Handle Out of Bounds Indices

Answer to the first post

I guess that is pretty much dependent on the problem that you are dealing with. Boundaries are always a problem that needs extra care.

In most applications, it would be an option to set those values to zero (and handle normalization factors of the filter appropriately).

Other options would be to reflect the data, so that the index -1 becomes 1, -2 becomes 2, and so on. Alternatively, there can also be cases where it would make sense to assume that it is periodic, so the index -1 would become N, -2 would bekome N-1 etc.

Hints regarding your update and follow-up doubt

You are right, that does not look right. Honestly, I havn't had a closer look at the formula you provided. However, I wonder if they are correct. I doubt that, honestly. The sums over $m$ and $n$, respectively, are just going in a positive direction. Having said that, it becomes apparent that you won't get a proper result, eg for $C^\prime(1,1)$, since the indices of $I$ are always below 0. That does not look right. Instead of coming up with lots and lots of formula here, I would like to point you to two links:

Quite visual explaination of 2D convolution/filtering

The obligatory Wikipedia link ;)

Nice visualization of the process, yet with bad audio

I like to think about the process in a rather visual way, as described in the first link. In my mind, I place the middle point of the filter ($K$) on the pixel I would like to evaluate. The other filter entries are then over the adjacent image pixels. Everything outside the original filter dimension is zero. Then you multiply this "extended" filter pixelwise with the image, and sum all up. That is the filter value at the given $(x, y)$ coordinate. You then move the filter's mid point by one pixel position and do the same thing again. Repeat this process until you have covered the full input image. See link #3


In a periodic square wave, how to change in a particular condition?

If $T_1$ remains constant and $T\rightarrow\infty$ then you get a single square pulse of length $2T_1$. There won't be another square pulse because the period is infinite.

$T_1$ constant and $T\rightarrow 0$ is indeed nonsensical because the definition of the signal assumes $T_1 < T/2$.

If $T$ is constant, you can get a constant function if $T_1=T/2$ (and not $2/T$).

If $T$ is constant and $T_1\rightarrow 0$ your function becomes 0. You will only get an impulse if you keep the area under the rectangle constant, i.e. if its height increases as its width decreases.


Techniques Used To Process Accelerometer Recordings

One suggestion (non-simple):  Obtain thousands of recordings from a broad spectrum of thousands of people engaged in your chosen activities to measure (with the activity externally logged, witnessed log book, video, etc.).  Feed all that data to several contemporary machine learning algorithms (DNN, etc.), and hope for a resulting inference method that meets your targets.


How can it be that there is a series of integrals in Fourier Series if it’s a projection on a continuous basis?

The short answer is that the sinusoidal components are only orthogonal if they are integer multiples of each other. Any other choice of frequency would not result in $<s1, s2> = 0$, which can be easily confirmed by comparing such cases (integer and non-integer related sinusoids). Thus the choice of basis signals is discrete in units of their frequency (the frequency domain), but the time domain waveforms that each of these discrete frequency values represent is continuous.

The frequency domain signal can still be a continuous signal if we allow zero-value signals, as it is just the non-zero values that would necessarily be discrete. That said, we can use discrete signals to describe the waveform in frequency without any loss of information.

The intuitive view starts from understanding that the Fourier Series Expansion states that we can represent a continuous time single-valued function (with no discontinuities) of duration $T$ as a series of sinusoidal components (sines and cosines) each with a frequency of $1/T$ and integer multiples.

Consider such an arbitrary function below:

Mathematically we need not limit time to $T$ in the the reconstruction; time can extend to $±\infty$ as long as we allow the signal in time to be periodic over period T, we will achieve the same result as derived in the Fourier Series Expansion. Since the waveform repeats exactly in time over every interval $T$, we need only express the result over period $T$ but extending it as such provides further insight into the Fourier Series Expansion, such as the inability to converge if there was a discontinuity between the end and the beginning of the interval. Here in the graphic I purposely made that boundary continuous so that the repeated waveform has no discontinuity at each transition, thus this is an example that the expansion will converge everywhere from $0$ to $T$ where in general such a continuity would not exist and you will typically see severe ringing at the start and end of the reconstructed waveform.

Thus the waveform components in the decomposition MUST also be cyclical over the interval $T$, otherwise they cannot possibly converge to the same waveform over the next period of $T$, meaning the signal cannot be periodic over $T$.  (As demonstrated with the possible first two components of the illustrated waveform as shown below.)

It is periodicity specifically in the time domain that is associated with discrete values of frequency in the frequency domain since the other frequencies if they did exist would destroy the periodicity over that interval: In order to be discrete in frequency, periodicity must exist in the time domain. Consider a sine-wave itself which is periodic it time and exists at one frequency. As demonstrated with the expansion an arbitrary periodic waveform can have multiple frequency components but they all must be integer multiples of the fundamental component at frequency $1/T$. I like to explain it like this: If I could call out your name exactly once per second, and exactly the same way each time--the resulting frequencies of my voice would only exist at 1 Hz and integer multiples of 1 Hz, all the way out to the 4KHz or so where spectral content exists.

This also works in reverse: periodicity in frequency is associated with discrete values in time (sampling). When we sample a time domain signal, we end up with multiple copies of the signal spectrum in frequency, spaced by the sampling rate. Same thing- in this case the frequency spectrum is continuous and periodic while the time domain signal is discrete. The DFT is an example that is discrete both in time and frequency, and therefore there is an implied periodicity both in time and in frequency!


While this is not much of an answer, here's a few comments while you wait for an expert to reply.

If this is one of your first DSP projects, it might be a bit ambitious just to dive straight into the fast implementations - you might want to start out with a naive implementation and then work your way up to the faster algorithms.

For the FFT-based approaches, you could check out how scipy does it:
http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.fftpack.dct.html

There's a lot of papers on the subject, but they might be a bit tricky to chew through. A couple of examples through some brief searches:

http://www.cin.ufpe.br/~vak/tg/papers/Chen-dct.pdf (1977 paper on a fast DCT)
http://arxiv.org/pdf/1203.3442.pdf (fast DCT-II)

It would probably be easier to understand from a book - I recall reading about the DCT and some fast algorithms with pretty elaborate explanations in "JPEG Still image data compression standard" by W. B. Pennebaker and J. L. Mitchell (I'm sure you can find loads of books on the DCT, specifically, though).

Also, this seems to be how they do it over at FFmpeg: http://www.ffmpeg.org/doxygen/2.1/dct_8c_source.html


Histogram equalisation (on its own) is not likely to work here because it corrects the brightness levels of the image as a whole.

What you are dealing with here is a uniform reduction in brightness akin to the effect of an optical filter.

What is really required is the correction of the mean brightness level over the "patch" that appears darker. Therefore, the problem now is to find the patch and correct it.

I am going to use Octave for the rest of this post but as you will see too it is not too difficult to transpose the idea to Python.

Let's have a look at the damage:

%Load the image and plot its average values by column.
a = imread('SomeImage.jpg');
imshow(a); %This should simply depict the image
%Let's take the cross section of the image by "squashing" all of its columns into their average.
q = mean(a);
plot(q);

This comes out like:

The most striking feature in this image is the sharp drop of the brightness level, about a third of the image width from the left.

Let's now look at the vertical cross section by "squashing" all of the image's rows into a mean column. All you have to add is:

q2=mean(a');
plot(q2);

And this comes out like:

In fact, if you have MATLAB somewhere handy, you might want to try to surf(a) and that will show you the "basin(s)" very clearly. For some reason, Octave's surf is getting its triangles tangled when I try to do it.

Right, so, what do we observe from this?

The image on either side of the discontinuities is predictable. It looks as if you simply took part of the waveform and slide it downwards...similar to the effect of a neutral density filter. There's simply a dip in the brightness that has a rectangular shape.
The problem on the first image is worse than the problem on the second image. On the second image we observe a sharp discontinuity which is approximately 4 pixels wide but it is very very sharp. On the first image, the discontinuity is about 8 pixels wide. Whatever data there was in those pixels has been "swallowed" by the gradient.
Working in 1 dimension (at first) you can see that just by either adding a value to the dip or subtracting a value from the "step", you can more or less match the means of the two regions on either side of the discontinuity. But for those pixels in the transition band, you are going to have to do something else (e.g. interpolation) to remedy the problem. The data in that region has been lost.

Where is the (big) patch? It starts approximately at 85,78 and extends all the way to the margins of the image. The X point is found by squashing all columns into one row and then looking at where the discontinuity starts and the Y point is found by doing the same but this time squashing all rows into one column.

Let's try to correct the X direction. First of all, we are going to get rid of the transition gradient. This is between pixels 85-93 (not across the whole image though, remember that this was a patch and that here we are trying to correct for the one dimensional case).

Q = q;
Q(85:93) = 0;
plot(Q);

Now, correct for the brightness. The difference is approximately 36. All I am doing here is to add a "step" (a constant value) from pixel 93 onwards, like so:

m=93;
Q2 = Q+36*[zeros(1,m),ones(1,length(Q)-m)];
plot(Q2);

This now looks like this:

What can you do for that gap of 8 values? You can fill it in with a constant value but that is going to show up as a "band" of 8 pixels and it will definitely be visible. The other thing you can do is to fill it with some Gaussian noise that has approximately pixel value 36 mean and enough variance so that you cannot distinguish it from the surrounding pixels. Unfortunately, because in this case the noise would be uncorrelated that would show up as banding again but this time much more "fuzzy". And finally, the other thing you could do is to train a model that would "generate" realistic values from the surrounding areas so that you "heal" the bands gracefully.

This is the basic idea. To automatically find the gradients you can either extract the gradient of the image with imgradient or do a similar calculation in a loop with a threshold on the local slope.

Your problems are not over here. As you can see in all images, there are more than 1 "basins" in the image. Your first image has a light band on the top right part which goes right into the darkness right below it.

Hopefully they are all rectangular so you can simply repeat the exact same process as many times as it is required until the image contains "no basins".

Hope this helps.


Different "zoom" after reconstruction from Radon transform using Fourier Slice Theorem

MathWorks' image processing toolbox describes the Radon transform as the projection of the image intensity along a radial line oriented at a specific angle. Now, let the set of parallel lines spaced 1 px apart is drawn in a fixed direction (let it be vertical lines) and you rotate the image bounded by the square 128x128 units (in which the phantom slice is inscribed) and for each 1 degree of rotation take the projections of the image (let it be horizontal projection of rotated image), integrating the image intensity (optical density, luminosity, etc) along the parallel lines of your fixed line set. Your object under test, in all its rotation orientations, is bounded not by the 128x128 square, but by a circle of 128·sqrt(2) diameter. To make sure that no region of the source image is missing, the Radon transform is performed for the circle of 128·sqrt(2) diameter.
Now, you perform the reconstruction of all the Radon transform data. The meaningful data are located in the circle of 128·sqrt(2) diameter, not in the 128x128 square, and you inscribe the reconstructed image into in the 128x128 square as is seen in the lines
x = -N/2:1:N/2-1; [X, Y] = meshgrid(x, x);
fhatcart = griddata(pX, pY, ghatpol, X, Y);
frec = ifftshift(ifft2(fhatcart));

of your code.
As the result, the structures of the source image are zoomed out by ~1.414.
You can adjust the meshgrid
x = -N/2:1.41421:N/2-1;

and the other parameters accordingly, or zoom in the reconstructed image in other ways, if you need to keep a 100% zoom.


How can I implement a triangular filter in MatLab, given it's Z-transform representation?

This MATLAB/OCTAVE code makes a causal implemtation the system $G(z)$ you mention. Note that $G(z)$ is non-causal, so you have to assume delays to make it causal. You can convert it to non-causal if you have the full data available before processing.

k = 5;
N = 2*k+1;

K = (k+1)^(-2);

b = K*[1, zeros(1,k) , -2, zeros(1,k), 1];
a = [1,-2,1];

[Hk,w]=freqz(b,a,1024*8,'whole');
figure,plot(w/pi,20*log10(abs(Hk)));title('frequency response');
ylabel('Magnitude (dB)');
xlabel('Frequency ( x \pi) ');

hn = impz(b,a,N);
figure,stem(0:N-1,hn);title('(causal) impulse response h[n]');


How do optical anti-aliasing filters work from a frequency domain perspective?

In Karel Fliegel, Modeling and Measurement of Image Sensor Characteristics, Radioengineering, vol. 13, no. 4, December 2004, he gives the optical transfer functions (OTFs) for different detector photosensitive area shapes. OTF is the Fourier transform of the spatial domain impulse response. The spatial domain impulse response, also called the point spread function (PSF), is a 2-d function that is constant-valued (assuming uniform sensitivity) inside the detector photosensitive area and zero outside it.

If there is an optical low-pass (blurring) filter in front of the sensor or implemented by shaking the sensor, another PSF describes it, and the total PSF is the convolution of the PSF of the sensor element and the PSF of the optical low-pass filter. The total OTF is the product of the two OTFs. Birefringenge-based optical filters (for example patent US 6937283 B1) create shifted copies of the image so their PSFs are sums of shifted and possibly scaled Dirac pulses. For example, a square arrangement of four Dirac pulses is equivalent to a cascade of a horizontal and a vertical comb filter.

An optical anti-aliasing filter also attenuates in-band frequencies that would not alias. Using a software sharpening filter the amplitudes of those frequencies can be restored. This also amplifies aliased frequencies, but because their attenuation in the combined OTF of the anti-aliasing filter and the detector is stronger, they retain some of that attenuation.

Here is an example analysis of a horizontal-vertical-separable total PSF that includes 1) a square photosensitive area, 2) a birefringent optical low-pass filter that has a PSF that consists of a square arrangement of Dirac pulses separated horizontally and vertically by 0.4284281154 times pixel width, and 3) a digital post-processing filter (operating on the pixel data), which can be separated to a cascade of a horizontal and a vertical filter, both with a PSF [-0.01536945896, 1.030738917, -0.01536945896]. A square grid of pixels with no gaps between pixels (100 % coverage of the imaging plane by photosensitive areas) is assumed. The analysis is done only along one of the separable dimensions for simplicity, but it would be better to do a 2-d analysis.

Figure 1. Red: combined one-dimensional PSF of the optical anti-aliasing filter, the square photosensitive area, and the digital sharpening filter. Turquoise: the ideal sensor PSF for a square grid (sinc function). Black: PSF of the square photosensitive area. Horizontal axis: displacement in units of one pixel width.

Figure 2. One-dimensional OTFs (Fourier transforms) of PSFs of Fig. 1, using the same color code. Horizontal axis: spatial frequency in units of radians per pixel width. There are some plotting artifacts that make parts of the red curve invisible.

Considering sinc (Fig. 1) as the ideal one-dimensional PSF, the mean square error for imaging spatial white noise is reduced by 40 % by using the optical anti-aliasing filter and the digital sharpening filter. The sum of squares error can be calculated as the integral over the square of the difference between sinc and the actual PSF, or by Parseval's theorem equivalently in the frequency domain using the differences of the OTFs (Fig. 2). I optimized the coefficients and the amount of separation in the optical anti-aliasing filter so that the mean square error is minimized.

Without the optical anti-aliasing filter, a digital post-processing filter (not a sharpening filter this time) with optimized coefficients [0.06765983302, 0.8646803339, 0.06765983302] only reduces the mean-square error by 11 % compared to a reduction of 40 % from the combination of an optical anti-aliasing filter and a post-processing sharpening filter, so it is apparent that optical anti-aliasing is useful.


Calculating the PDF of a waveform from its samples

What I'd go with is essentially Jason R's "random resampler", which in turn is a presampled-signal based implementation of yoda's stochastic sampling.

I've used simple cubic interpolation to one random point between each two samples. For a primitive synth sound (decaying from a saturated non-bandlimited square-like signal +even harmonics to a sine) it looks like this:

Let's compare it to a higher-sampled version,

and the weird one with the same samplerate but no interpolation.

Notable artifact of this method is the overshoot in the square-like domain, but this is actually what the PDF of the sinc-filtered signal (as I said, my signal is not bandlimited) would also look like and represents the perceived loudness much better than the peaks, if this were an audio signal.

Code (Haskell):

cubInterpolate vll vl v vr vrr vrrr x
    = v*lSpline x + vr*rSpline x
      + ((vr-vl) - (vrr-vll)/4)*ldSpline x
      + ((vrr-v) - (vrrr-vl)/4)*rdSpline x
     where lSpline x = rSpline (1-x)
           rSpline x = x*x * (3-2*x)
           ldSpline x = x * (1 + x*(x-2))
           rdSpline x = -ldSpline (1-x)

                   --  rand list   IN samples  OUT samples
stochasticAntiAlias :: [Double] -> [Double] -> [Double]
stochasticAntiAlias rs (lsll:lsl:lsc:lsr:lsrr:[]) = []
stochasticAntiAlias (r:rLst) (lsll:lsl:lsc:lsr:lsrr:lsrrr:t)
    = ( cubInterpolate lsll lsl lsc lsr lsrr lsrrr r )
          : stochasticAntiAlias rLst (lsll:lsl:lsc:lsr:lsrr:lsrrr:t)

rand list is a list of random variables in range [0,1].


Average spectrum representation - SonicVisualiser

It would not be useful to average in time because vocalisations are represented directly by their spectrogram. (For example: http://www.u.arizona.edu/~ohalad/Phonetics/notes/Formants%20Spectrograms%20and%20Vowels.PDF or even http://www2.coloradocollege.edu/dept/ev/Research/Faculty/OVALItems/FordOsborneVocabulary/_SouthernVocabularyTable.html)

I would suggest that you export your spectrogram as an image and then isolate and layout the areas you are interested in, in some other software like GIMP. For more information, please see: http://www.sonicvisualiser.org/doc/reference/1.9/en/#imageexport and http://www.gimp.org/

Hope this helps.


Plotting the magnitude response of a biquad filter

Someone over at the KRV forum pointed out that I need to normalize my coefficients so that a0 = 1.

b0 /= a0;
b1 /= a0;
b2 /= a0;
a1 /= a0;
a2 /= a0;
a0 = 1;

Solved my problem.


How channel state information is calculated from Sounding Packet

According to pages 5-7 in this, the Channel State Information (CSI) contains 2 measurements per matrix element.

The CSI essentially tells the Access Point (AP) two things: where is the client (bearing) and how far away on that bearing it is. However, because of multipath propagation, the "bearing" might not correspond to the true bearing of the client but the direction at which the signal is "louder". This means that the signal could be bouncing on a wall and then finding the client, rather than traveling to it in a straight line.

So, the task we have now is to calculate these two numbers. We assume that the channel identification pseudorandom sequences in the pre-amble are transmitted at the same time from the AP:

To estimate the $A$, simply sum the "strength" of the received signal at each antenna.
To estimate the $T$, you need to estimate the time differences between the times the signal was received. At this point, you need to decide on the ordering of the elements.

You can estimate both using a matched filter whose $h$ is the pseudorandom sequence.

The matched filter will return a relatively high output when it detects the pseudorandom sequence at its input.

So, imagine a bank of matched filters behind each receiver. Each one of the antennas receives the AP signal at slightly different times and amplitudes. Each one of the matched filters produces a slightly different amplitude (the $A$ component) when it "sees" the pseudorandom sequence at its input. Consequently, because of the difference in arrival times, each one of these high outputs will occur at slightly different times. To estimate $T$, take one of your antennas as a reference element and measure the time that passes between two high outputs at each matched filter.

You might also find reading up on the fundamental of the following techniques:

The Matched Filter and its applications in Radar and Sonar
Beamforming

Which sooner or later will lead you to Adaptive Beamforming which is what you are dealing with here.

Phased Arrays

Which, again, is what you are dealing with here.

Pseudorandom binary sequences (PRBS)

Which is what your header sequence is right here

PRB Sequences with applications to channel characterisation

In brief, by examining the distortions that a PRBS undergoes when transmitted through a "channel" (whatever this might be), you get an estimate of the channel's condition or state which you can then use to adapt some aspect of transmission (Forward error control sequence, beamforming, etc).
Characterisation of an acoustic communication channel using PRBS
Transmitter identification using PRBS
And more...

Hope this helps.


Spectral estimation in Spectrum analyzer

The methods you're mentioning are methods for estimating spectra from a (typically discrete) time-domain signal – what the "classical" spectrum analyzer (SA) does is not that!
The classical sweeping SA consists of four parts:

A tuneable oscillator
a mixer
a good filter
a power detector

Here, 1. and 2. are used to mix down different parts of the spectrum - either by sweeping (typical for completely analog SA) or stepping through a set of frequencies.
The output of the sweeping mixer is then filtered by 3., which effectively means 1.-3. are a method of getting a varying piece of RF passband down to a frequency where 4. can be used to measure how much power is in that piece of RF spectrum.
You'd use the instantaneous frequency of the oscillator as x-coordinate, and the power "seen" by the power detector as y-coordinate. That's your point!
Note how you're not looking at the whole spectrum at once. You'll find "probability of intercept" considerations in literature, i.e. how likely is it that a bursty signal is actually active when the SA "tunes by" and observes the frequency they're on.
There's nothing here that really deals with the time-domain signal. For all that matters, the signal simply "has" some energy at some frequency. No need to estimate the spectrum from some time signal – they're actually measuring the energy at a single frequency, then move on to the next.
If you want to know more, there's really many introductions to using a SA, and you'll find the same components (oscillator, mixer, filter, power detector) again, and how you can configure them (sweep speed/dwell time, amplification, filter bandwidths, sensitivities).
Now, I didn't use "classical" with quotation marks without reason: For the last ~20 years, SA typically do have the possibility to sample the mixed-down signal and analyze it to improve the resolution, to implement better filters, or to do signal analysis beyond estimating a spectrum. Typically, these measurement methods are pretty limited, even in upper-class SA; they are literally labeled "FFT" without specifying exactly what they do, but you can expect things to be Bartlett's method for as long as the SA dwells on one frequency. How the SA pieces together different spectral estimates from different tuning frequencies is something that I haven't even found in manuals. Armed with an SDR and time you can figure that one out, but in all honesty, that's rarely what you really need an SA for: SDRs and even oscilloscopes can give you a better spectral estimates (though I can tell you, as DSP guy, the nomenclature of some oscilloscope vendors describing their spectral estimation methods drives me nuts).


The quadratic surface is determined by the autocorrelation matrix of the data, which is always positive definite or positive semi-definite. This means that any stationary point is always a minimum. In the worst case, this minimum is not unique if the matrix is singular, but it can never be a saddle point.


FM synthesis algorithm is not very efficient

I think you're overthinking the integrator:

Since your $s$ is a discrete signal in itself (you're doing DSP!!), the integral $S(T) = \int\limits_0^T s(t) \,\mathrm dt$ of it at any time $T$ actually collapses to the running sum!

So, as you've noticed, when you've got the value for any time $T$, the value at time $T+T_s$ (one sample later) is simply $S(T+T_s) = S(T) + s(T+T_s)$.

This should really be the least complex step in your FM.

Now, when you think about real-world FM radios, you'll find that the simple circuits (but also the more complex ones) don't actually contain an actual integrator (which you could build with an Operational Amplifier and a capacitor, for example).

What they really contain, and what makes the spectral side of things much nicer to understand, is a low-pass filter (you'll notice an integrator acts as low-pass filter). You'll notice there will be other filters in real-world broadcasting systems, so-called preemphasis filters. I'd call them a historical artifact from when receivers couldn't be trusted to equalize overly well, and that leads to a noise shaping that leads to decreased SNR (much noise) at higher audio frequencies, if you don't pre-amplify these selectively.

Now, the classical FM transmitter thus consists of an integrator-alike low-pass filter, and a preemphasis filter with a  slight high-pass characteristic. You could probably combine these two.

My conclusion to that topic always is that FM is super easy to explain superficially, but looking deep, it gets really ugly, and the way the sausageaudio is actually produced at the receiver is loosely defined by the imperfections of 70 years old radio electronics and the technological abilities of that time. On an analytical level, FM gets especially cumbersome to deal with; try yourself on the derivation of the Bessel-shaped spectra of a single audio tone in FM.  This, in fact, lead to us (as in: the communications theory basic course at university) dropping FM from the curriculum alltogether: if you'd build an audio transceiver system in this time and age, you can always get better and more resilient than FM utilizing the same power·bandwidth product by going digital. Thus, we try to avoid weighing down students with hard math and questionable applicability, even if it means that they aren't taught how the most prolific broadcasting standard of the world works.

This really shouldn't stop you, though! This is a worthwhile experiment, and building an FM transmitter is actually pretty satisfying soon as it works. Note that real-world receivers are pretty resilient to choosing wrong filters; we've actually had a bug in GNU Radio for what I guess was more than a decade in which the filter coefficients for the de-emphasis (ie. receiver side) filter were plain wrong, and no-one noticed. People simply don't expect exactly great-sounding audio over FM.


How to simplify an exponential functions in MATLAB?

Symbolic math is included in MATLAB, but rarely used in signal processing. There are cases, nevertheless, when we wish to use symbolic math packages to simplify some expressions. One such alternative is also Wolfram-Alpha web application.
Anyway, apart from a misunderstanding, you indeed have a (simple) mathematical problem:
Given a complex number $z$ , you have
$$ z + z^* = 2 ~\mathcal{Re}\{ z \} $$
or a complex-valued function $f(t)$
$$ f(t) + f^*(t) = 2 ~\mathcal{Re}\{ f(t) \} $$
If you look at your function $g(t)$, you will see that it's a sum of a complex-valued function $f(t)$ and its conjugate $f^*(t)$, where
$$f(t) = 14.9 ~j~e^{j 31^o} ~ e^{j 2t} $$
and its conjugate
$$f^*(t) = - 14.9 ~j~e^{-j 31^o} ~ e^{-j 2t} $$
hence
$$g(t) = f(t)+f^*(t) = 2 \mathcal{Re}\{ f(t) \} $$
Then you can see that
$$\mathcal{Re}\{ f(t) \} = \mathcal{Re}\{ 14.9 ~j~e^{j 31^o} ~ e^{j 2t} \} =\mathcal{Re}\{ 14.9 ~e^{j (121^o +2t)} \} = 14.9 \cos(121^o +2t)$$
Then your $g(t)$ is:
$$ g(t) = 2 \times \mathcal{Re}\{f(t)\} = 2 \times 14.9 \cos(121^o +2t) $$


What is the precise definition of a nonlinear signal?

In my opinion, the term nonlinear signal doesn't make much sense. If you refer to nonlinear signal processing then this is not about the processing of nonlinear signals, but about the nonlinear processing of signals, i.e., about the use of nonlinear systems.

A system $\mathcal{T}$ is nonlinear if it doesn't obey the superposition principle. I.e., if $y_1(t)=\mathcal{T}\{x_1(t)\}$ is the response to input $x_1(t)$, and $y_2(t)=\mathcal{T}\{x_2(t)\}$ is the response to input $x_2(t)$ then

$$\mathcal{T}\big\{\alpha x_1(t)+\beta x_2(t)\big\}=\alpha\mathcal{T}\big\{x_1(t)\big\}+\beta\mathcal{T}\big\{x_2(t)\big\}$$

is generally not satisfied.

Well-known examples of such nonlinear systems are median filters and Volterra filters, but there are many trivial operations, such as taking the magnitude of a signal, which must be classified as nonlinear processing.


Area under waveform in t-domain including 2 delta functions

please someone confirm if following  devlopment is right or  wrong .becuase answer of exercise question isn't provided in book "F F kuo"

we know

$f(t)\xrightarrow{\mathcal{F}}F(j\omega); $ where, $F(j\omega) =\displaystyle\int_{-\infty}^{\infty} f(t) \ e^{-j\omega t} \ dt\implies F(0)=\text{area under time domain curve }$

So, $(a)\ F(j\omega)=2+T.\ Sa{(\omega T/2)}+Ke^{-j\omega T}\implies F(0)=2+T+K=\displaystyle\int_{-\infty}^{\infty}f(t)dt$

but area is given zero therefore , $K=-(2+T)$ .

and it is resonable also becuase strength of impulse at $t=T$ must be equal to the negative summation of  area of regions preceding it.

similarly,

for part $(b)\ $   $K=-T$ because time is starting from $t=0_{+}$ that's why excluding impulse of strength  $2$ at $t=0$

note :

$Sa (\omega T/2)=\dfrac{\sin\left(\dfrac{\omega T}{2}\right)}{\omega T/2}$


The error should be minimized between the equalized signal and what? There is no known sequence in this case.

But there is: if the channel is still good enough so you can decode a packet successfully, you can reconstruct the original signal by simply re-encoding the packet.


Delay signal according to normal distribution

Normal distribution regardless of whatever mean and variance you provide, can still produce large positive and negative numbers(albeit with a small probability).

Do you want to specifically use normal distribution, or all you need is a random delay generator? If latter is the case, you can use rand() from MATLAB, also you could use randi if you need integer delays.

Also you could use something like this to generate uniformly distributed random numbers in a specific range.


Method of reconstructing a band-limited signal from discrete samples

Some recent cell phone models use something like a Cirrus Logic CS42xx series audio IO chip, which seems to use a digital polyphase interpolation filter, a sigma delta modulator, followed by a switched capacitor DAC and low-pass filter.

Sinc interpolation (or, given finite hardware, a polyphase FIR kernel similar to a windowed Sinc) is one high quality method of reconstruction from bandlimited samples.

But reconstruction isn't the main issue with respect to audio quality.  Quality or fidelity is more likely constrained by the audio data compression algorithm used (or required by standardization), and how much spectrum bandwidth is allocated by the comunications standard and regulations.


Is it feasible to try and classify an enclosed space as being a. Very crowded b. Less crowded & c. Sparsely crowded using image processing?

This sounds feasible (and fun!), but it will be hard to get very accurate.  My guess is that, like your tags suggest, this would be a combination of edge detection / space detection to get a measure of density and face/hand detection to get a measure of human content.  You could then train up an SVM or RandomForest model on the data, tweak, and repeat.


Audio: How to measure dynamic range (in sense of dynamic compression/ "loudness war")

One important of them is the dynamic range (in sense of "loudness war"). By which in this context is not meant the SNR but the (naive said) "variance" of the "occuring energy".

Loudness War is neither dynamic range nor SNR nor variance of occuring engery (although it's related). It's mostly about peak limiting.
For better or for worse, there is a strong desire in the industry to make sure that your song is percieved louder when being played back to back with somone else's song on the radio, a streaming service or a record company executive review meeting.
The song is delivered as a fixed point audio file and this format has a maximum sample amplitude that it can represent. To make it "loud" you need to shove as much energy as possible into the file without any one sample exceeding the the allowable max amplitude. In practice this is done the following way

You record and mix in an unconstrained format (floating point).
The last step is mastering, where the mix session is turned into a fixed point file with whatever format is required (sample rate, bit depth, labelling, streaming constraint, etc.)
During mastering so-called "peak limiting" is applied. You raise the overal level of the song and then shave off the few samples that exceed the maximum allowable amplitude. This generates distortion.
Repeat this until you have achieved the maximum "tolerable" distortion and you can't make it louder without adding too much crap.

In essence you create a trade off between overall loudness and "acceptable" distortion. "Loudness war" basically means that this trade off has skewed hard to towards "loud" and "lots of distortion". A good example to analyze is Tom Petty's album "Hypnotic Eye" the CD version is peak limited, but the high res version is not, so you can actually compare directly what damage the peak limiter has done.
This is quite painful to listen to especially since it's mostly unneccessary. The original version is not as loud, but I can simply adjust the volume to compensate. However, the distortion in the peak limited version is always there.
Back to the original question: you want to look at the "RMS to Peak", since that's exactly what the peak limiters is trying to reduce.  Since these days the "peak" in a file is almost always very close to unity (unless it's "mastered for Itunes"), you can simply take the RMS of the file. Here are some typical examples

0 dB: Square wave
-3 dB: Sine wave
-7.8 dB: "Maggie's Farm" by Rage Against the Machine
-10.8 dB: Averabe of the Billboard 2013 top 100
-15 dB: Average of 10000+ "typical" songs including a wide mix of genre's and release years

Many streaming services don't like this variability since it requires the customer to keep "riding" the volume control and therefore prescribe mastering levels (example -14 dB LUFS which is based on an ITU norm). They may remaster the content or ask the labels to do it which the labels may or may not do. As a result, the same song when played from Spotify may be mastered differently than the same one played from Apple Music.


I am not sure that you would call this a practical application as much as a practical implication, but knowing (or forcing) a signal to be odd or even enables you to process it in a certain way.

For example, if you know (or force) a signal to be even, then the imaginary part of its Discrete Fourier Transform (DFT) goes to zero and you can represent this signal just as a sum of $\cos(\cdot)$ functions. This relates the DFT to the Discrete Cosine Transform (DCT).

Specifically:

$$X[k] = \sum_{n=0}^{N-1} x[n] \cdot e^{-i \frac{2 \pi k n}{N}}$$

This is a compact way of saying:

$$X[k] = \sum_{n=0}^{N-1} x[n] \cdot \left ( cos\left ( \frac{2 \pi k n}{N} \right ) - i \sin\left ( \frac{2 \pi k n}{N} \right ) \right ) $$

This makes it clear that $X \in \mathbb{C}$ (omplex) and is composed of a real part (where $x[n]$ is multiplied by the $\cos(\cdot)$) and an imaginary part (where $x[n]$ is multiplied by the $\sin(\cdot)$).

But, $\cos(\cdot), \sin(\cdot)$ are even and odd functions respectively and the way their symmetry works when combined with integration means that if a signal is an even periodic function, then it can be represented by a series of cosines.

OK, but I only have a sound recording that could be anything...how do I know (or force) it to be an even function? (or an odd function).

One of the key assumptions of the Discrete Fourier Transform is that $x[n]$ is a periodic function. That is, even if you pass some random recording of a person counting "one two three four", all that the DFT sees at its input is "one two three four one two three four one two three four one two three four...". In other words, once you reach the length of $x[n]$ (which is $N$), you wrap back to its start.

So, how do we make this an even signal?

Notice that in even signals $f(x) = f(-x)$. We could therefore generate a new signal that is double the length of $x[n]$ and once we reach the end of $x[n]$ we repeat it by playing it backwards. In this way, $x[N+n] = x[N-n]$ or "one two three fourouf eerht owt eno one two three fourouf eerht owt eno one two three fourouf eerht owt eno..." and so on.

In terms of a "demo", the DFT of [1,2,3,4] has an imaginary component while the DFT of [1,2,3,4,3,2] doesn't.

I do not think that you will find a direct application as such, it will probably be something that "exploits" the implication that a function is even or odd.

Hope this helps.


What is the decay rate (time constant) of two One-Pole LPF's cascaded in series?

If both filters have the same time constant, i.e., if their individual impulse responses are
$$h(t)=e^{-\alpha t}u(t),\qquad \alpha>0\tag{1}$$
then the impulse response of the cascade is
$$h_{tot}(t)=te^{-\alpha t}u(t)\tag{2}$$
If both filter have different time constants with individual impulse responses given by
$$h_i(t)=e^{-\alpha_it}u(t),\qquad \alpha_i>0\tag{3}$$
the total impulse response is given by
$$h_{tot}(t)=\frac{e^{-\alpha_1t}-e^{-\alpha_2t}}{\alpha_2-\alpha_1}u(t)\tag{4}$$
I don't think it makes much sense to ask for an "effective" time constant of the cascaded systems. How should one define such a time constant? Of course, for large $t$ the impulse response $(2)$ decays as $e^{-\alpha t}$, and the impulse response $(4)$ decays as $e^{-\min(\alpha_1,\alpha_2)t}$, but that doesn't really say much about the systems' behavior.


Approximate the magnitude response of an analog filter with a digital filter, starting from its gain expression

Did you try to plot this curve in matlab.
I am not sure whether it is any standard high,low or band pass filter.
So, to calculate the impulse response there might be two ways.
1. take the inversefft of the NdB and then do the convolution between your impulse input and the inversefft signal of filter usinc conv(x,y) command.
2. Or convert the impulse input into frequency domain by fft transform in matlab and then just multiply the gain and fft of input impulse and then convert them in time domain.
Can you do it using above method?
Do you want a matlab code for this method?


Do I get less accurate results if I do not use a power of 2 points for my FFT?

Assuming a proper implementation, and discarding the small rounding differences that you'll get just by doing a different sequence of operations, no, there's no substantive drop in accuracy if you slightly pad your DFT length to the next power of 2. Whether it ends up actually being any faster is dependent upon your implementation, however. Why is that?

It's a very common misconception that there is a single Fast Fourier Transform algorithm and that you must use a power-of-2 size for best performance. If you look for a description of an FFT algorithm, you'll often see the radix-2 decimation-in-time or decimation-in-frequency techniques explained, likely because they're the easiest to illustrate. However, even the seminal FFT technique, the Cooley-Tukey algorithm, generically factorizes the FFT size into smaller numbers, not just powers of 2.

Using a good FFT library, you'll get the best performance if your FFT size can be factored into a number of small prime factors. The FFT library will then have optimized implementations of DFT kernels for each of these primes, which can then be recombined appropriately to yield the full set of DFT outputs. As I alluded to in a comment on another question, modern libraries will often give good performance for all prime factors ~13 and below.

With that said, you may find that by padding up to the next power of 2, it's possible that your transform might become slower. If you're already using an FFT size that your library implementation is well-suited for, you're just adding extra work for yourself by padding the size out. The best way to judge this is to benchmark a few candidate sizes and see which does best on your platform. If you have to make an automated choice of a good FFT size, then based on the characteristics of the radixes that your library supports, you can choose the next size that has an appropriate set of prime factors.


Converting Audacity Filter Curve EQ into transfer function and applying it to a signal via python

You can use the manage/export functionality to convert your Audacity filter preset into a text file. The file contains frequency/value pairs plus some side info. It looks like Audacity builds an 8192 FFT grid and uses a B-spline to interpolate the frequency/value pairs over the entire grid. I assume that just do an inverse FFT plus some circular time shifting to build a linear phase impulse response and than do something like fftfilt() to apply the filter.
You can just read the text file into python, parse it, and then follow the same steps.


I guess that would be a slew rate limiter. This is concept is mostly used an amplifier design as a practical constraint of the circuit. I haven't seen it applied as a digital filter. It is certainly very non-linear and it would make a poor "smoothing filter" as it's highly dependent  on the absolute amplitude. Could you shed some light on the specific application and explain what a regular low pass filter can't do the desired smoothing job?


Using a Kalman filter to predict samples between measurements

I'm sorry that it was not very clear to me what you were asking. But I think I have developed some sense on it as follows:

You are physically capturing 10 noisy sensor measurements per second, which I call it as the physical measurements. But for some reasons, you want to fill in between those true-physical data samples with new synthetic data, generated mathematically at a rate of 100 samples per second (10 times more dense than available with the physical sensor measurement)

The solution to this problem has different names, based on your interpretation and which method you choose to follow, such as a sampler, a  predictor, an estimator, a curve fitting or an interpolation...

It seems, you want to use a Kalman filter (an estimator actually) to generate those new synthetic samples in between the true physical measurements. And your question is whether you should feed those synthetically generated new sampes back into the Kalman filter (for residual calculations) as if the data were generated from a true physical sensor, so that the filter can use it in residual calculations to get rid of the noise present in it and refine the state estimates?

That's what I understand from your question set up and my answer to it is no you cannot. Those new data cannot replace true phsyical measurements. They do not involve new information that's called as the innovations process of a random process (that you're trying to estimate), hence they cannot be used to make refined computations or decisions alone.

In fact, the Kalman filter will only propagate the current state $X(t)$ into the next phsyical measurement point $X(t+T_s)$ by using the system dynamic matrix alone. (i.e., it will solve the differential equation numerically based on the initial condition provided from the previous state estimate, and sample this solution at the requested high data rate)


DFT and integer valued basis functions'sf requencies

Could you develop on the usage of integer frequencies?

When dealing with sampled discrete signals, there can only be an integer amount of resolvable frequencies within a captured interval.

When you take an $x(t), t \in \mathbb{R}$ and apply sample-and-hold on it at some sampling frequency $Fs$, it is turned into an $x \left[n \cdot \frac{1}{Fs} \right ], n \in \mathbb{N}$.

Theoretically, when $t \in \mathbb{R}$, you can "interrogate" any time instance or interval of $x$. But after sampling, when $n \in \mathbb{N}$, the smallest interval you can interrogate from your signal is $\frac{1}{Fs}$ or between successive samples $n, n+1$.

Even if you tried to sample faster after the sample-and-hold, the value you would get would be the last known value from the last time the sampling took place.

Therefore, within a block of $N$ samples, you get $N$ resolvable frequencies. Even if you tried to evaluate the Discrete Fourier Transform in more than the available $N$ frequencies, all that you would get would be interpolated values between the $N$ resolvable bins of it.

Which brings us to what Stanley Pawlukiewicz remarks: The fact that $n$ is integer, does not mean that you can only represent integer frequencies. Only that you can represent a fixed amount of them.

Just as you can only "see" $n$ and then $n+1$ discrete time instances of $x$, so you can only resolve some $f, f+ \frac{Fs}{N}$ discrete frequencies within it, where $N$ is the amount of samples you have collected.

  ...possibly relate it to the DCT where it seems to me that, e.g. the second basis vector is like half a cycle (as shown on my plot), and half a cycle is like a fractional frequency?

I hope that it is clear from the above discussion that "...half a cycle..." can well be N samples of some $f$ sampled at the right $Fs$ (?).

Hope this helps.


How to shift a signal on the $y$-axis to match another signal in MATLAB

Matlab's polyfit function can be used to do linear regression. This can be used to find a linear relationship (offset and gain) that maps one dataset to the other.


Block normalization of histogram of oriented gradient

I stucked at this question too, because overlapped blocks will be divided twice, right?

and so I found these :

1 http://www.geocities.ws/talh_davidc/#cst_extract

2 https://stackoverflow.com/questions/32417531/hog-what-is-done-in-the-contrast-normalization-step

to verify my answer, I done this:
on 1, they using image size HOG126x63(), block size is 18x18, with 50% overlapping blocks means it had 9px stride.
126/9, 63/9 = (13, 6) blocks.

they mentioned that final vector would be 6318 = 78*81 == 78*9*9 (78*number_of_cells_in_one_block*number_of_bins)
turns out, 13*6 = 78, so final vector is
(blocks_in_image_col *blocks_in_image_row *number_of_cells_in_one_block *number_of_bins)

conclusion is, one block is counted as one set vectors (2*2*9 as in image) , moving the blocks will produce one set of vectors again. there won't be cells normalized twice (it will make the histogram spike at first and flat on remaining due to recursive dividings)


How to normalize the power of a QAM signal?

In $16$-QAM, four constellation points
$\left(\pm\frac{d}{2},\pm\frac{d}{2}\right)$ are "interior" points while
the remaining $12$ are "edge" points or "corner" points for which
at least one of the coordinates is $\pm\frac{3d}{2}$.
Thus, in
$n$ measurements $(X_i,Y_i)$ of successive transmitted symbols,
some of the $X_i$ will equal $\pm\frac{d}{2} + NI_i$ while
others will be $\pm\frac{3d}{2} + NI_i$.  Similarly,
some of the $Y_i$ will equal $\pm\frac{d}{2} + NQ_i$ while
others will be $\pm\frac{3d}{2} + NQ_i$.  Here, $NI_i$ and $NQ_i$
are independent zero-mean Gaussian random variables with variance
$\sigma^2$ representing the noise. Note that we can expect
about $\frac{3}{4}$-th of the
samples to have at least one large coordinate. Now, consider
that if we want to estimate the value of $d$ from the measurements,
and want to use the average energy calculations, then
$\sum_i (X_i^2+Y_i^2)$ is one possibility, but we are weighting
equally the symbols with large means and the symbols with small
means, and the noise component is a more significant part of
the symbols with small means than it is of symbols with large
means.  It might be better to use a weighted sum to estimate
$d$ (or $d^2$), and one extreme possibility is to give zero weight to
interior points ($(X_i$ and $Y_i$ are both small) and to use
only the large values of $X_i$ and $Y_i$ (edge and corner
points) in the estimation.  If the number of samples is large,
this could be as simple as taking those (approximately) $3n/4$
samples in which either $X_i$ or $Y_i$ (possibly both) is
large, and averaging the absolute values of the largest $X_i$ and
$Y_i$ to get an estimate of $3d/2$.

Once upon a time, when
analog circuits were in use, these kinds of calculations would
have been difficult to implement in hardware, and this effect
carries over into some present-day implementations that faithfully
follow the methods of the past that were originally devised
because of hardware limitations.
But these days, with digital computation being ubiquitous,
perhaps something new can be attempted.


Essentially, you are evaluating the partial sum of a Fourier Series of the rectangular periodic pulse. The pulse is symmetric around $t=0$, so the sine-parts can be omitted, leaving only the cosine-parts, as you have shown. However, you miss an essential thing: the frequency zero (i.e. the DC part), which is the overall mean of your waveform.

I dont know exactly, how you came up with the different values for amp, but it turns out the waveform with dcyc=0.5 has an amplitude of approximated 0.78. So, the mean of a waveform with a given $d_{cyc}$ is given by

$$a_0=0.78d_{cyc}+(1-d_{cyc})(-0.78)$$

You need to add this mean value to your Fourier Series to generate the equal-amplitude waveform. See below:

length = 1024
nharms = 8
dcyc = 0.2
step = 2*np.pi/length
harm = 1

A = 0.78
mean = (dcyc) * A + (1-dcyc)*(-A)

j = np.arange(4*length)

table = 0
for i in range(nharms):
    amp = np.sin(harm*dcyc*np.pi)/harm
    table = table + amp*np.cos(step*harm*j)
    harm = harm + 1
table = table + mean

plt.plot(j, table);
plt.grid(True)


The proper way to define your frequency vector after a DFT is as follows. Let $N$ be your DFT length, and $f_s$ be your sampling rate in Hz. Furthermore, define an $N$-length frequency vector $\bf{f}$, where each element $f_i = i$, for $i = 0, 1, 2, ... N-1$.

Now your frequency vector in hertz is simply going to be $\bf{f}$$\frac{f_s}{N}$

Now, assuming you are DFT'ing a real sequence, simply pick all elements with frequency values less than $\frac{f_s}{2}$, and you are in business!

You can also see my answer here for actual code.


Can I compute Fourier series without looping through all frequencies?

Make sure you do not mix up FFT and DFT, they gives identical results, but are different.
DFT complexity is about $O(N^2)$ while FFT complexity is $\frac{34}{9}N\log_2(N)$. This makes a huge difference is computation time.

I do not know any faster algorithm for computing a complete spectrum than the FFT. If you need a subset of the spectrum (only few frequency bins), you could look at Goertzel Algorithm.

According to wikipedia, as a rule of thumb, Gortzel algorithm will be faster than FFT if :
$$ M < \frac{5N_2}{6N}\log_2(N_2)$$

Where

$M$ : Number of frequency bin you need
$N$ : Number of point in your dataset
$N_2$ : The first power of 2 bigger than N.


This is always going to require a lot of computation, especially if you want to process as many as 2000 points.  I'm sure there are highly-optimized solutions for this kind of pattern-matching already, but you have to figure out what it's called in order to find them.

Since you're talking about a point cloud (sparse data) instead of an image, my cross-correlation method doesn't really apply (and would be even worse computationally).  Something like RANSAC probably finds a match quickly, but I don't know much about it.

My attempt at a solution:

Assumptions:

You want to find the best match, not just a loose or "probably correct" match
Match will have some small amount of error due to noise in measurement or computation
Source points are coplanar
All source points must exist in target (= any unmatched point is a mismatch for the entire profile)

So you should be able to take a lot of shortcuts by disqualifying things and decrease computation time.  In short:

pick three points from source
search through target points, finding sets of 3 points with the same shape
when a match of 3 points is found, check all the other points in the plane they define to see if they're a close match
if more than one match of all points is found, choose the one with the smallest sum of 3D distances error

More detailed:

pick a point from the source for testing s1 = (x1, y1)
Find nearest point in source s2 = (x2, y2)
d12 = (x1-x2)^2 + (y1-y2)^2
Find second nearest point in source s3 = (x3, y3)
d13 = (x1-x3)^2 + (y1-y3)^2
d23 = (x2-x3)^2 + (y2-y3)^2

for all (x,y,z) test points t1 in target:
    # imagine s1 and t1 are coincident
    for all other points t2 in target:
        if distance from test point > d12:
            break out of loop and try another t2 point
        if distance ≈ d12:
            # imagine source is now rotated so that s1 and s2 are collinear with t1 and t2
            for all other points t3 in target:
                if distance from t1 > d13 or from t2 > d23:
                    break and try another t3
                if distance from t1 ≈ d13 and from t2 ≈ d23:
                    # Now you've found matching triangles in source and target
                    # align source so that s1, s2, s3 are coplanar with t1, t2, t3
                    project all source points onto this target plane
                    for all other points in source:
                        find nearest point in target
                        measure distance from source point to target point
                        if it's not within a threshold:
                            break and try a new t3
                        else:
                            sum errors of all matched points for this configuration (defined by t1, t2, t3)

Whichever configuration has the least squared error for all other points is the best match

Since we're working with 3 nearest neighbor test points, matching target points can be simplified by checking if they are within some radius.  If searching for a radius of 1 from (0, 0), for instance, we can disqualify (2, 0) based on x1 - x2, without calculating the actual Euclidean distance, to speed it up a bit.  This assumes that subtraction is faster than multiplication. There are optimized searches based on more arbitrary fixed radius, too.

function is_closer_than(x1, y1, z1, x2, y2, z2, distance):
    if abs(x1 - x2) or abs(y1 - y2) or abs(z1 - z2) > distance:
        return False
    return (x1 - x2)^2 + (y1 - y2)^2 + (z1 - z2)^2 > distance^2 # sqrt is slow

3D Euclidean distance is $d = \sqrt{(x_1 - x_2)^2 + (y_1 - y_2)^2+(z_1 - z_2)^2}$
but square roots are slow. Since we're just comparing relative distance, I think we could drop the sqrt and consider "distance" to be sums of squares.

The minimum calculation time would be if no 2-point matches are found.  If there are 2000 points in the target, this would be 2000 * 2000 distance calculations, though many would be disqualified by a subtraction, and the results of previous calculations could be stored so you only have to do $2000 \choose 2$ = 1,999,000.

Actually, since you will need to calculate all of these anyway, whether you find matches or not, and since you only care about nearest neighbors for this step, if you have the memory it's probably better to pre-compute these values using an optimized algorithm.  Something like a Delaunay or Pitteway triangulation, where every point in the target is connected to its nearest neighbors.  Store those in a table, then look them up for each point when trying to fit the source triangle to one of the target triangles.

There's a lot of calculations involved, but it should be relatively fast since it's only operating on the data, which is sparse, instead of multiplying lots of meaningless zeros together like cross-correlation of volumetric data would involve.  This same idea would work for the 2D case if you found the centers of the dots first and stored them as a set of coordinates.


What is the difference between the PSD of a deterministic and stochastic signal?

Generally, one cannot compute the Fourier transform of any random process, because of integrability reasons. The PSD is one way to define it, based on autocorrelation.

And now, this mathematically sound definition of a periodogram is consistent in some way: if you have a realization of a random process, on a limited time interval, then its Fourier transform can be computed, and is consistent with the PSD, which you'd expect.


Symbol Timing recovery for modulation producing ISI

In general the ISI is never severe enough for the symbol timing recovery algorithm not to work, but it does degrade its achievable performance on the $P_{be}$ vs. $\dfrac{E_b}{N_0}$ curve.  In other words, for the same error rate, you need better SNR to make up for the ISI.

ISI will degrade the S-curve of the Timing Error detector, causing it to flatten a bit and reduce the TED gain.  Gaussian noise will do the same.


Why we need Frequency and Time information at the same time

For the source, go to end of the answer

Suppose one day you got one note which has some thing written to it, say "Major frequency components are 10 Hz, 25Hz, 50 Hz and 100 Hz". Somehow, you understood that its time-series representation is a very important thing (may be master-piece work of a great musician, or some national security matter, anything). So you decided to find it.

You start by plotting the spectrum which looks like below:

(There can be some noises/disturbances in between these frequencies, but most people neglect it)

Now you want to find its time-series representation. Since we have worked it a lot in our signal processing class in college, most probably, we get below signal:

$$
x(t)=cos(2 \pi 10 t)+cos(2\pi 25 t)+cos(2 \pi 50 t)+cos(2 \pi 100 t)
$$

Great !!!

But is this the only time series representation for those frequencies?

Consider another time series signal which plays a sine wave of frequency 100 Hz for 300 ms, then 50 Hz sine wave for next 300 ms, then 25 Hz sine wave for 200 ms, then 10 Hz sine wave for 200 ms. A total of 1000 ms long signal. If you plot it, it will look like below:

Actually above signal also have similar frequency spectrum (shown below):

Those ripples are because of sudden changes in frequency. Amplitude changes are because of change in how long each frequency was played. They are not of importance to us right now.

So what to do now? How will you pick one? Even if we change the order they were played, still we get the same frequency spectrum. So how many such combinations you have?

  In short words, without knowing when those frequencies were played, it
  is very difficult to find its time series representation, and that is
  where, wavelets come into picture. It provides not only the
  frequencies present in the signal, but also, the time it was present
  in the signal.

UPDATE (17-Feb-14)

Question : How do you relate that we require both spatial (time) and frequency information in an image at a same time ?

Answer : Me too, had the same question before, so after seeing your comment, I thought I should find one example (I haven't studied wavelets in the context of images before wavelets). And below explanation is the summary of an answer in Quora.com (all credits to that author, Ron Reiter) which I felt as a good example.

You might be knowing that we use the DCT for jpeg compression of images. It preserves the low frequency data while throws away some of the high frequency data in the images, based on the image quality we need.Or in other words, this method is good with low frequency data, but not that accurate with high frequency data.

For example, in an image with sky and a tree, the areas of sky will be preserved accurately while some details of tree will be thrown off.

In DCT, we take 8x8 blocks of images and encode it for all frequency of data. For very low frequencies, we could have used higher block size while for very high frequencies, we may have to use lower block size than this.

Now apply our wavelet method as explained in 1D case above. What if we can find both frequency and their location in the image? We can encode the low frequency areas in big chunks while high frequency areas in low chunks. So it produces better results.

Please check the link for more details

SOURCE:

THE WAVELET TUTORIAL by ROBI POLIKAR

All the images, explanation, everything, are extracted from this link. It is a must-read introduction of wavelets to beginners

Amara's Wavelet Page
Wavelet in image compression - Ron Reiter at Quora.com


Do we use closed loop or open loop information in Bode plot, Nyquist plot and Root Locus?

The bode plot is just a plot showing the frequency representation of your system. Any transfer function has one, hence you can use it to see the response of both closed or open loop. Nyquist and the root locus are mainly used to see the properties of the closed loop system. The root locus shows the position of the poles of the c.l. system as the gain of your controller changes. While nyquist diagram contains the same information of the bode plot. But instead of splotting two real graphs of the amplitude and the phase of the complex function the nyquist plot shows the same function in the complex plane


Nevermind, I figured it out. For my purposes, all I really want is the cepstrum of the energy-normalized signal, which after working through the algebra is equivalent to adding a delta function to the original cepstrum. (See below.)

For a real, discrete, finite-time signal $\vec{s}$, we can compute its energy as $E(\vec{s})=\vec{s}\cdot\vec{s}$. To energy-normalize the signal (i.e. ensure that $\text{energy}=1$), we multiply $\vec{s}$ by the scalar $E(\vec{s})^{-1/2}$.

The cepstrum of the the signal $\vec{s}$ is defined as $\text{cep}(\vec{s}) = \mathfrak{F}^{-1}\big(\text{log}\big(\big|\mathfrak{F}(\vec{s})\big|^2\big)\big)$.
Therefore the cepstrum of the energy-normalized signal can be computed as: $\text{cep}\big(E(\vec{s})^{-1/2}\cdot\vec{s}\big)\\ = \mathfrak{F}^{-1}\big(\text{log}\big(\big|\mathfrak{F}\big(E(\vec{s})^{-1/2}\cdot\vec{s}\big)\big|^2\big)\big)\\ =
\mathfrak{F}^{-1}\big(\text{log}\big(\big|E(\vec{s})^{-1/2}\cdot\mathfrak{F}(\vec{s})\big|^2\big)\big)\\ =
\mathfrak{F}^{-1}\big(\text{log}\big(E(\vec{s})^{-1}\cdot\big|\mathfrak{F}(\vec{s})\big|^2\big)\big)\\ =
\mathfrak{F}^{-1}\big(\text{log}\big(E(\vec{s})^{-1}\big)\cdot\mathbf{1} + \text{log}\big(\big|\mathfrak{F}(\vec{s})\big|^2\big)\big)\\ =
\mathfrak{F}^{-1}\big(\text{log}\big(\big|\mathfrak{F}(\vec{s})\big|^2\big) - \text{log}\big(E(\vec{s})\big)\cdot\mathbf{1}\big)\\ =
\mathfrak{F}^{-1}\big(\text{log}\big(\big|\mathfrak{F}(\vec{s})\big|^2\big)\big) - \mathfrak{F}^{-1}\big(\text{log}\big(E(\vec{s})\big)\cdot\mathbf{1}\big)\\ =
\text{cep}(\vec{s}) - \big[\text{log}\big(E(\vec{s})\big)\;\; 0\;\; 0\;\; \dots\;\; 0\big]$

Where $\mathbf{1}$ denotes the vector of all 1's.


What effect does rotation in the spatial domain has on phase in Fourier transforms?

This question has a simpler answer for the 2-D continuous-space Fourier transform but itsdDiscrete Fourier transform based verification requires some elaboration and careful implementation as @MarcusMüller has already mentioned.

On the continuous-space it can be shown that rotating an image $f(x,y)$ by $\theta$ radians CCW (counter-clock wise) on the $xy$ plane, will also rotate the corresponding CSFT (continuous-space Fourier transform) $F( \Omega_1 , \Omega_2 )$ by $\theta$ radians CCW :

CCW rotation by $\theta$ on the $xy$ plane can be described by the following transform:

$$
x' = x \cos(\theta) - y \sin(\theta) \\
y' = x \sin(\theta) + y \cos(\theta) \tag{1}
$$

where $x,y$ are the original coordinates and $x',y'$ are the rotated (new) coordinates. A corresponding inverse rotation by $\theta$, (or CCW by $-\theta$ ) has the following transform :

$$
x = + x' \cos(\theta) + y' \sin(\theta)  \tag{2} \\
y = - x' \sin(\theta) + y' \cos(\theta)
$$

Given the original image $f(x,y)$, we rotate it CCW by $\theta$ to obtain the rotated image $g_r(x,y)$, then the following will hold :

$$ g_r(x,y) = f( x \cos(\theta) + y \sin(\theta) , - x \sin(\theta) + y \cos(\theta) ) \tag{3} $$

Note that in Eq. (3), the arguments of the original function $f(\cdot,\cdot)$ are described by the inverse rotation relation to coordinates $x,y$ of $g_r(\cdot,\cdot)$.

The CSFT of the rotated image $g_r(x,y)$ is:

$$
\boxed{ G_r(\Omega_1,\Omega_2) = \iint_{-\infty}^{\infty} g_r(x,y) e^{-j( \Omega_1 x + \Omega_2 y) } dx dy} \\
$$

$$
\begin{align}
&= \iint f( x \cos(\theta) + y \sin(\theta) , - x \sin(\theta) + y \cos(\theta) ) e^{-j( \Omega_1 x + \Omega_2 y) }   \tag{4} \\
&= \iint f( x' , y' ) e^{-j[ \Omega_1 (x' \cos(\theta) - y'\sin(\theta)) + \Omega_2 (x' \sin(\theta) + y' \cos(\theta)) ]} dx' dy'  \tag{5}\\
&= \iint f( x' , y' ) e^{-j[ ( \Omega_1 \cos(\theta) + \Omega_2  \sin(\theta) ) x' + (-\Omega_1 \sin(\theta) + \Omega_2 \cos(\theta)) y' ]} dx' dy'  \tag{6}\\
&= F(  \Omega_1 \cos(\theta) + \Omega_2  \sin(\theta) ,   -\Omega_1 \sin(\theta) + \Omega_2 \cos(\theta) ) \tag{7}
\end{align}
$$

In moving from step (4) to (5), make the substitution defined in Eqs (1) and (2), and note that the Jacobian of the transformation is unity. See page 778 of Calculus Adams 6E for an explanation and derivation of this result.

In moving from step (6) to (7) we simply recognise the expression in Eq.(6) as the CSFT of the original signal $f(x,y)$ with arguments given according to relation in Eq.(3) given as below:

$$G_r(\Omega_1,\Omega_2) = F(  \Omega_1 \cos(\theta) + \Omega_2  \sin(\theta) ,   -\Omega_1 \sin(\theta) + \Omega_2 \cos(\theta) ) \tag{8} $$

Eq.(8) is analogous to Eq.(3) and states that the CSFT of the rotated image is also roated by the same amount;  $\theta$ CCW. This, therefore, means that phase is also rotated by the same amount.

Verification of this result within the discrete-space sequences $f[n,m]$ and their corresponding DFT sequences requires specific care given to potential aliasing caused by inadequate sampling of the rotated image $g_r(x,y)$, and  proper circular extension of the rotated sequence $g_r[n,m]$.

Without going into mathematical details, I'll put the results of the following implementation of the rotation using Matlab/ Octave (rotation is wrt upper left corner)

We can still see artefacts of non-ideal implementation...

Note that with typical images it's harder to realize the ideal result.


Don't understand the relation about the wireless communication flow chart and TDMA illustration diagram

TDMA wraps around the outside of your block diagram of the communication link.  There are $M$ users who are buffering the data that they wish to transmit. The TDMA switch at the input selects each user in turn and connects its buffer to the block marked Source Coding for a fixed amount uf time. During this time, the contents of the user's buffer are source coded, channel coded, transmitted etc.  When the user's time is up, the TDMA switch connects the next user's buffer to the Source Coding block, etc.  At the receiver end, there is a synchronous switch that routes the Source Decoder's output to the appropriate user's buffer which the user empties at a lower rate (if it so chooses so that it gets its data at a constant rate rather than in bursts).


Measured signal from tachometer has ripples, why and how to minimize?

Well, regarding the why: don't know the underlying mechanism of your tachometer; there's many, many ways to build one.

Minimization: compared to your signal of interest, these seem to be rather high in frequency. So, a reasonably designed low-pass filter will help. Don't forget to account for group delay, should you need to feed back the resulting signal in e.g. a control system.


The plot that shows the sequence $x[n]$ is also called as the lollipop plot and it shows the individual samples of $x[n]$ with their weights and locations.

Each lollipop represents a discrete-time unit-impulse signal, $\delta[n]$, which is located at the position it's shown and weighted by the value of the signal there. The superposition sum of all those impulses compose the signal $x[n]$.

Now, your signal has 3 lollipops located at $n = \{-1,0,1\}$ and weighted by $\{1,2,3\}$ respectively.

Since an impulse located at the position $n=d$, and weighted by $a$, is given by
$$ a \cdot \delta[n-d]$$

by using this, you can compose the signal $x[n]$ from its lollipops as:

$$ x[n] = 1 \cdot \delta[n-(-1)] + 2 \cdot \delta[n-(0)] + 3 \cdot \delta[n-(1)] $$

$$ x[n] = 1 \cdot \delta[n+1] + 2 \cdot \delta[n] + 3 \cdot \delta[n-1] $$

Note that a generalization of this also yields the sifting-sum decomposition of any signal $x[n]$ as

$$x[n] = \sum_{k=-\infty}^{\infty} x[k] \delta[n-k] $$


Convolution - Massive spike in noise section

If your searched-for signal is partially similar to itself at some offset, that can also show up when correlated against an actual input signal containing it.

Test your searched-for signal by autocorrelation against itself to check for this problem.  If you want a sharp convolution peak against test input, you need a searched-for signal that has very low self-similarity at any lag.


Is MATLAB documentation wrong on when to use discrete versus continuous wavelet transform?

[EDIT] Let us start from a verbose and simple example. In $\mathbb{R}^2$, the Mercedes-Benz pattern constitute a tight frame: redundant, energy-preserving:

Each vector $x=[x_1,x_2]$ yields three three projections $[p1,p2,p3]$, from which you can recover $x$ in a number of robust ways: if the projection value are corrupted, rounded, even if one is cancelled, you can still recover $x$, at least approximately. Because of the redundancy $2\to 3$ in the representation. This is not the case with the orthogonal basis the $\{[1,0],[0,1]\}$: if you lose one projection, $x$ is almost lost (you only keep one direction).

In other words, you don't need all coefficients with a frame of vectors (or wavelets). This basic example extends to more generic frames or bases. You can find more on that within the topic of oversampled filter banks, and the paper Filter bank frame expansions with erasures. Those are used for redundant channel coding, for instance.

In the continuous setting (both the signal and the wavelet transform), due to the huge redundancy, many possible inverses exist. Notably, given an analysis wavelet, the requirements on the synthesis wavelets to verify the admissibility conditions (on a pair of wavelets) are very mild. On these topics, I would suggest two references:

Matthias Holschneider, Wavelets, an analysis tool, 1995
Bruno Torrésani, Analyse continue par ondelettes, 1995

The topic is quite involved, may involve a lot of complicated functions spaces, different for the wavelets and the signals. A wavelet analysis is given by:

\begin{equation}
W_f(b,a)\ =\ \int_{\mathbb{R}}\ f(t)\ \tfrac{1}{\sqrt{a}}\,\psi^{*}(\tfrac{t-b}{a})\ d
t\ .
\end{equation}
Provided that $\psi$ is admissible, i.e. when the
two constants $c^{\pm}_{\psi} =
2\pi \int_{0}^{+\infty}\tfrac{|\hat\psi(\pm\omega)|^{2}}{\omega} d\omega
< \infty$ are finite and equal\footnote{When $\psi$ is sufficiently
  regular, this condition reduces to a zero-average requirement, that is,
  $\int_{\mathbb{R}} \psi(t)\ d t = 0$}, that is,
$c^+_\psi=c^-_\psi=c_\psi < \infty$, the signal $f$ may be recovered
from the coefficients $W_f(b,a)$:
\begin{equation}
f(t)\ =\
\tfrac{1}{c_{\psi}}\int_{0}^{+\infty}\!\!\!\int_{\mathbb{R}}\ W_f(b,a)\
\psi_{(b,a)}(t)\ d b\,\tfrac{d
a}{a^{2}}.
\end{equation}

You can find about the same idea in Torrésani's book, page 20:

Still considering a continuous signal with finite energy, but trying to build a discretized  wavelet frame, you can sample in the wavelet space with perfect reconstruction provided the sampling is dense enough with respect to some lower bound computed from the wavelet. So effectively, you don't need all (continuous) values. But the number of samples needed could be denombrable.

Basically, you have a frame when for every unit energy signal $f$ projected onto the set of frame elements $\psi_n$:

$$A\le \sum |<f,\psi_n>| \le B\,.$$

If $A=B=1$, you recover the case of orthogonal expansions. Ingrid Daubechies (for one) showed, from the continuous formula above, that one can find some discretizations $a_n$ and $b_n$ to build frame elements

$$\psi_n = \psi(\tfrac{t-b_n}{a_n})$$

and recover $f$ from the $|<f,\psi_n>| $, sometimes with iterative formulae that converge faster if $A$ is close to $B$.

Moving on to discrete signals (where differentiability is not well-defined), you can build unions of wavelet basis, or wavelet frames like the stationary wavelet transform. In this redundant (but with a finite number of samples)  expansion, you can pick a subset of samples, that corresponds to, at minimum, a specific choice of basis vectors. Once again, you do not need all samples.

THe rationale is that if an $N$ sample signal get represented by $M> N$ coefficients, when they are well formed (frame coefficients for instance), different inverses exist, some of then discarding parts of coefficients.

For a basis (non redundant) expansion however, all samples are required in general.


How to adjust receiver gains to avoid saturation and quantization noise to optimise post digital processing?

To avoid quantization or limit cycle oscillations, you can add hysteresis to your AGC decision loop (a form of "debouncing").


okay, for step 1, let's get the continuous time math down right and simplified to as much as possible.  i am changing some notational convention slightly to be more consistent with Electrical Engineering practice:

$$v(t) = \sum_{k=0}^{N-1} V_\mathrm{step}(1-e^{-(t-k\Delta t_\mathrm{step})/\tau})u(t-k\Delta t_\mathrm{step})$$

$$V(j\omega) =\mathscr{F}\Big\{v(t)\Big\} = V_\mathrm{step}\left(\sum_{k=0}^{N-1}e^{-j\omega k\Delta t_\mathrm{step}}\right)\left(\frac{1}{j\omega} +\pi\delta(\omega) - \frac{1}{\frac{1}{\tau}+j\omega}\right)$$

there is a well known closed form for the summation of a geometric series:

$$ \sum_{k=0}^{N-1} x^k = \frac{x^N -1}{x-1} $$

so

$$\begin{align}
 \sum_{k=0}^{N-1} e^{-j\omega k\Delta t_\mathrm{step}} &= \sum_{k=0}^{N-1} (e^{-j\omega \Delta t_\mathrm{step}})^k \\ \\
 &= \frac{(e^{-j\omega \Delta t_\mathrm{step}})^N -1}{e^{-j\omega \Delta t_\mathrm{step}}-1} \\ \\
 &= \frac{e^{-j\omega N \Delta t_\mathrm{step}}-1}{e^{-j\omega \Delta t_\mathrm{step}}-1} \\ \\
 &= \frac{e^{-j (\omega/2) N \Delta t_\mathrm{step}}(e^{-j (\omega/2) N \Delta t_\mathrm{step}}-e^{j (\omega/2) N \Delta t_\mathrm{step}})}{e^{-j(\omega/2) \Delta t_\mathrm{step}}(e^{-j(\omega/2) \Delta t_\mathrm{step}}-e^{j(\omega/2) \Delta t_\mathrm{step}})} \\ \\
 &= e^{-j (\omega/2) (N-1) \Delta t_\mathrm{step}}\frac{(e^{j (\omega/2) N \Delta t_\mathrm{step}}-e^{-j (\omega/2) N \Delta t_\mathrm{step}})/(2j)}{(e^{j(\omega/2) \Delta t_\mathrm{step}}-e^{-j(\omega/2) \Delta t_\mathrm{step}})/(2j)} \\ \\
 &= e^{-j (\omega/2) (N-1) \Delta t_\mathrm{step}}\frac{\sin(\omega N \Delta t_\mathrm{step}/2)}{\sin(\omega \Delta t_\mathrm{step}/2)} \\ \\
\end{align}$$

so i think your final expression in continuous frequency of the spectrum should be

$$V(j\omega) = V_\mathrm{step}\frac{\sin(\omega N \Delta t_\mathrm{step}/2)}{\sin(\omega \Delta t_\mathrm{step}/2)} \left(\frac{1}{j\omega} +\pi\delta(\omega) - \frac{\tau}{1+j\omega \tau}\right) e^{-j (\omega/2)(N-1)\Delta t_\mathrm{step}} $$

does that look right to you?


how gabor filter is used on image texture edge detection

A periodic geometrical texture may appear as a localized wiggling pattern, with dominant orientations and 2D frequencies. Its correlation with Gabor atoms is bigger when the atom has orientation and frequencies close to the texture's, thus yielding feature vectors allowing texture detection or classification.


How do Anti-Aliasing Filters filter the signal?

Yes that is correct. The width of the transition band is inversely proportional to the length (or memory) of the filter which means to say that you would need an infinite amount of time to achieve your perfect filter. Therefore for practical reasons you decide how much aliasing would be tolerable (similar in many ways to deciding how many decimal places you need to include for a number describing the length of something) and then design your filter accordingly. The tighter the rejection the longer will be the filter delay and higher will be the overall complexity.


Handling a number of simultaneous measurements in Extended Kalman Filter?

If your measurements are batch, so is your filter.

You can do a number of things.

You can sequence your measurements in natural order and update the filter for each.

You can do the previous and then approach the problem as a fixed lag smoother.

There is a literature on out of sequence and delayed updating as well, like:

  Ensemble Kalman filter for multisensor fusion with multistep delayed
  measurements Sirichai Pornsarayouth; Masaki Yamakita

  2011 Aerospace Conference Year: 2011 Pages: 1 - 10 Cited by: Papers
  (1) IEEE Conferences


LTI system, Fourier transform, unit step function, sine wave

I can see three problems here:

1) $x(t)$ is a cosine wave but only where $t$ is positive, not $\omega$

2) You must keep the step function when applying Euler's identity, and I can't see the cosine frequency in your second expression for $x(t)$ either. Always make sure your identity is right, with WolframAlpha for example

3) Do not mix up your $\omega$ that is a parameter to $\cos(\omega t)$ and the $\omega$ variable of the frequency domain in $X(j\omega)$. Use different notations, like $\cos(\omega_0t)$


How does kernel work in unsharp masking?

The first kernel that you gave,

$$
\left[ \begin{array}{ccc} -1 &&-1 &&-1 \\
-1 && 8 && -1 \\
-1 &&-1 &&-1 \end{array} \right]
$$

is interpreted as follows: for each pixel in the image, take 8 times its value (hence the 8 in the center of the kernel) and then subtract the values of all the 8 neighboring pixels (hence all of the surrounding -1 values). This is a form of a highpass filter, which tends to accentuate edges.

Think about what you get if you were to run this kernel over a patch of a constant color. If the region has a uniform color, then it doesn't contain any edges. The kernel above would yield a value of zero when applied to a 3-by-3 pixel block that has the same color, which is desirable. Not only is the filter highpass, it also completely attenuates zero frequency, which is also good.

If you change the center value to 14 instead of 8, it no longer has the property where it will output zero when applied to a region of constant color. Essentially, the kernel will have a different spatial frequency response, which leads to the differences in the second output image.


Least Angle Regression (LARS) without Matrix Inversion

If you want to solve for single value of $ \lambda $ in the model:

$$ \arg \min_{x} \frac{1}{2} {\left\| A x - b \right\|}_{2}^{2} + \lambda {\left\| x \right\|}_{1} $$

Then you can use Coordinate Descent method which is the fastest and simplest and doesn't require any matrix inversion.

I have a MATLAB code for in my $ {L}_{1} $ Regularized Least Squares Analysis GitHub Repository:


Difference between these two integrators

There are different methods to approximate integration in discrete time. The most straightforward ones are the forward and backward Euler methods, and the trapezoidal method. A discrete-time system with transfer function

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

implements the forward Euler method.

Similarly, the backward Euler method is implemented by

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

Finally, the trapezoidal method is implemented by

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

These discrete-time integrators are also described on this mathworks page.


Possibility of using Power allocation with Maximum likelihood estimation

What you mention could be right in low modulation order, such as QAM but in 16-QAM and 32-QAM, the amplitude is changed also, I mean you will have different amplitude levels, so you can't implement that detector.


Dependence of $\tt n$ in $\tt fft(data,n)$ in auto correlation of data

Doing an $N \times M$ point FFT for that sized image to calculate the autocorrelation will result in time-aliasing.

To multiply two signals of length $N$ in the frequency domain, the FFT length must be $2N-1$.


The problem comes from sampling the phase. For the given phase function, the changes in the phase from one value $t_i$ to the next one $t_{i+1}$ can become too large, so phase unwrapping will not work properly. If you simply removed the factor $2\pi$ in the exponent, the problem would be solved already (note that the factor shouldn't be there anyway according to definition of $y(t)$ in your first formula):

s = exp(1i*phi);

If you want to keep the factor $2\pi$, then use more sampling points.


Formulas of the Fourier transform family

As a complement to the above answers, we can also classify the transformations based on the orthogonality condition(s) used to derive each pair of transformations.

Fourier Transformation:
$$\frac{1}{2\pi}\int^{+\infty}_{-\infty}e^{ikx}dx = \delta(k)$$

Fourier Series and DTFT:
$$\frac{1}{2\pi}\int^{\pi}_{-\pi}e^{ikx}dx = \delta_{k,0}$$
$$\frac{1}{2\pi}\sum^{+\infty}_{k=-\infty}e^{ikx} = \delta(x),\;\;\;\;k\in Z, x\in (-\pi,\pi)$$
I group Fourier Series and DTFT together, because of the duality. The continuous, periodic function in time domain of Fourier Series, can be seen as the frequency domain function in DTFT. Therefore, there is not much difference between the two, at least on mathematical ground.

DFT:
$$\frac{1}{N}\sum^{N-1}_{n=0}e^{i\frac{2\pi}{N}kn} = \delta_{k,0},\;\;\;\;k=0,1,...N-1$$

Note: Take caution when dealing with the Delta functions.

Kronecker Delta: It can be treated as a discrete elementary function, and commonly used under a summation.
Dirac Delta: It cannot be derived by elementary function, although many introductory texts use $$\delta(x-a) = \lim_{n\to\infty}\sqrt{\frac{n}{\pi}}e^{-n(x-a)^2}$$ to explain it. And treat it as an ordinary function with a value of $\infty$ at $0$, and a value of $0$ otherwise. Delta function is a distribution and has it meaning definite when used under an integral.


why only impulse signal for convolution?

Because the impulse response completely characterizes an LTI system. The reason is an arbitrary input $x(t)$ can be written as an infinite sum of time shifted and scaled Delta functions:

$$x(t) = \int_{-\infty}^{+\infty}x(\tau)\delta(t-\tau)d\tau$$

or alternatively, using the convolution operator, $$x(t) = x(t)*\delta(t)$$

Properties of linearity and time invariance lead to the following sequence of input-output pairs:

$\delta(t) \longrightarrow h(t)$
Time invariance: $\delta(t-\tau) \longrightarrow h(t-\tau)$
Linearity: $x(\tau)\delta(t-\tau) \longrightarrow x(\tau)h(t-\tau)$
Linearity: $\int_{-\infty}^{+\infty}x(\tau)\delta(t-\tau)d\tau \longrightarrow \int_{-\infty}^{+\infty}x(\tau)h(t-\tau)d\tau$

The last input-output pair can be written as $$x(t) \longrightarrow y(t) = x(t)*h(t)$$

This way, any output can be computed using the input $x(t)$, the impulse response $h(t)$, and the convolution integral (or sum, in discrete time).

@a concerned citizen gave you some other hints regarding your question.

On the contrary, if you use the unit step function, you get

$$\int_{-\infty}^{+\infty}x(\tau)u(t-\tau)d\tau = \int_{-\infty}^t x(\tau)d\tau$$

or alternatively, using the convolution operator, $$\int_{-\infty}^t x(t) = x(t)*u(t)$$

Once again, just for comparison, you get the following sequence of input-output pairs:

$u(t) = \int_{-\infty}^t \delta(\tau)d\tau \longrightarrow s(t) = \int_{-\infty}^t h(\tau)d\tau$
Time invariance: $u(t-\tau) \longrightarrow s(t-\tau)$
Linearity: $x(\tau)u(t-\tau) \longrightarrow x(\tau)s(t-\tau)$
Linearity: $\int_{-\infty}^{+\infty}x(\tau)u(t-\tau)d\tau \longrightarrow \int_{-\infty}^{+\infty}x(\tau)s(t-\tau)d\tau$

where $s(t)$ is the step response, the system output when the unit step is presented at the input. The last input-output pair can be written as $$\int_{-\infty}^t x(\tau)d\tau \longrightarrow x(t)*s(t) = x(t)*u(t)*h(t)$$

The latter can be written as $$(x(t)*h(t))*u(t) = y(t)*u(t) = \int_{-\infty}^t y(\tau)d\tau$$

So after all, what you get at the output is the integral of the output you would get using the delta decomposition of an input signal. That is, your output must be differentiated to obtain the actual system output for a given $x(t)$.


Topological skeleton matching techniques

I think you should look/google in the direction of graph matching.

This is because you can represent your skeleton as a graph by creating nodes at endpoints and crossings, and creating edges for the connected endpoints and crossings. There is much more research done in the direction of graph matching than just topological skeleton matching. However for the latter I found this paper, which you might find interesting.

You can also look through the references of the Wikipedia page of Topological skeletonization, there are also some references about matching skeletons.