Nyquist with only poles at the origin I need to solve this three exercises:
$$G_1(s) = \frac{1}{s}$$

$$G_2(s) = \frac{1}{s^2}$$

$$G_3(s) = \frac{1}{s^3}$$

This is what I did:

For the single pole ($1/s$), magnitude of A is $\infty$ and phase is 0

For the single pole ($1/s$), magnitude of B is $\infty$ and phase is -45

For the single pole ($1/s$), magnitude of C is $\infty$ and phase is -90

For the 2 poles one ($1/s^2$), magnitude of A is $\infty$ and phase is 0

For the 2 poles one ($1/s^2$), magnitude of B is $\infty$ and phase is -90 ($=2*(-45)$)

For the 2 poles one ($1/s^2$), magnitude of C is $\infty$ and phase is -180 ($=2*(-90)$)

For the 3 poles one ($1/s^3$), magnitude of A is $\infty$ and phase is 0

For the 3 poles one ($1/s^3$), magnitude of B is $\infty$ and phase is -135 ($=3*(-45)$)

For the 3 poles one ($1/s^3$), magnitude of C is $\infty$ and phase is -270 ($=3*(-90)$)

MATLAB can't do these Nyquist analysis in a correct way. That's why I'm asking here.

Thanks

Your answer is correct. Each pole gives 90 degrees phase lag. Furthermore, here can you download a file lnyquist1 which is able to give a more accurate result when you have a pole at zero.

num = 1;
den = [1 0 0 0];
lnyquist1(num,den);

In addition, WolframAlpha is also able to give an accurate result, see here.


FIR Filter length Upper Bound Say I want to design a low pass filter with stop and pass band ripple constrained to $\epsilon$ and transition region $\delta$, what is the smallest FIR Filter length that suffices? Can you please provide a derivation?

Many people have asked themselves this question. Unfortunately, there is no exact solution, just approximations. One relevant paper is Practical Design Rules for Optimum Finite Impulse Response Low-Pass Digital Filters by Herrmann, Rabiner and Chan, where they basically fit a formula to their data obtained from lots of designs. A simplified formula is [1]

$$M=\frac{-10\log_{10}(\delta_1\delta_2)-13}{2.324\;\Delta\omega}$$

where $M$ is the filter order, and $\Delta\omega=\omega_s-\omega_p$ is the width of the transition band.

Also the Matlab function firpmord.m uses some formula to estimate the necessary filter order, but they also say that the required order is often underestimated.

[1] Discrete-Time Signal Processing, Oppenheim, Schafer, Buck (2nd ed.), p. 502

Try the FDAtool in MATLAB to experiment with various possible solutions.


Determine whether the system is a FIR or IIR by looking the transfer function I have the following system:
$$
y[n]=\frac{1}{3}(x[n+1]+x[n]+x[n-1])
$$
After the Z-Transform we get

$$
\frac{y[z]}{x[z]}=\frac{z^2+z+1}{3z}
$$
which is of course the transfer function of the system. Now, I know that it would be a FIR if $\sum h(z) < \infty$ with $h(z)$ the impulse response. But As far as I know, it should be clear if it is FIR or IIR by just looking the transfer function. How can I know that?

The forward $z$-transform is written with negative powers on $z$. So, $z$ transform of your difference equation would be:

$$
Y(z) = \frac{1}{3}\left(zX(z) + X(z) + z^{-1}X(z)\right)
$$

yielding a transfer function of:

$$
H(z) = \frac{Y(z)}{X(z)} = \frac{1}{3}(z + 1 + z^{-1})
$$

The easiest way to determine whether a filter is IIR or FIR is to identify its pole locations. For FIR filters, there is a rule for this that is based on the structure of the impulse response:

If the system is causal (i.e. it is zero for all $n < 0$), then it is FIR if all of its poles are located at the origin ($z=0$).
If the system is anticausal (i.e. it is zero for all $n \ge 0$), then it is FIR if all of its poles are located at $z=\infty$.

Your system is neither causal nor anticausal. Instead, it has a causal component ($\frac{1}{3}(x[n] + x[n-1])$) and an anticausal component ($\frac{1}{3}x[n+1]$). So, it can be helpful to look at each component separately:

$\frac{1}{3}(x[n] + x[n-1])$: The $z$-transform of the causal portion of the transfer function is $\frac{1}{3}(1 + z^{-1})$. This expression $\to \infty$ as $z \to 0$, so it has a pole at $z = 0$. This meets the criterion for an FIR response.
$\frac{1}{3}x[n+1]$: The $z$-transform of the anticausal portion of the transfer function is $z$. This expression $\to \infty$ as $z \to \infty$, so it has a pole at $z = \infty$. This meets the criterion for an FIR response.

This would imply that you gave a difference equation for a non-causal FIR filter.

Transfer function with only zeros is FIR.
If both zeros and poles are present, then it is IIR.
Please avoid pole zero cancellation for this identification method.


How to choose a windowing function when performing a Chirp $\mathcal Z$-transform? I've got some data sampled at 32 Hz. The only interesting stuff is happening between about 0-8 Hz and I've got to do real time batch processing and detect events with a delay of no more than a second. So that means I can only perform FFTs with 32 bins or less, to remedy the lack of resolution I'm using a Chirp $\mathcal Z$-transform (czt in Matlab).

I was wondering if there were special considerations to be made when you choose a windowing function on the original time data before when performing a CZT instead of an FFT?

When we add the window function to the definition of the DFT we get

$$ X(k) = \sum_{n=0}^{N-1} w[n]x[n] e^{-j \omega n} $$
with $\omega_k= \frac{2\pi k}{N}, k\in \mathbb{Z} , $ $w[n]$  the window function and $x[n]$ the data. The multiplication of the data $x[n]$ with the window function $w[n]$ in the time domain corresponds to a convolution in the frequency domain.

The CZT and the FFT are basically just different ways to compute this sum.  Each algorithm has its advantages and disadvantages. The FFT is fast, the CZT is more flexible. If you use both algorithms  to compute the same coefficients $X(k)$ you get the same results. So, in this case, the window choice considerations should be the same for CZT and FFT.

So far I can answer only have of the original question: If the CZT is used to compute the same frequency bins as delivered by the FFT, there is no difference with respect to window selection.

Note1: Are you aware of the issues regarding frequency resolution of different DFT algorithms, as discussed e.g. in this thread ? The frequency resolution topic is the same for "CZT vs. FFT" as for "zero padding of FFT" or for "Goertzel vs. FFT".

Note2: This was the basic message of my original post.

If  the CZT is used to compute the DFT on a finer frequency resolution, things are more complicated. In this thread the effect on the DFT is discussed. To make it short, the CZT interpolates the DFT. This interpolation makes window selection more involved.

If (as intended by the OP) the CZT is used to interpolate the DFT at a finer frequency resolution, the underlying maths remain the same, however the different frequency resolution must be kept in mind.

This is illustrated in a very simple example. The first figure shows the absolute values of a rectangular window. The FFT (without zero padding or similar things) computes the DFT right at frequencies where the window is zero, here marked in orange.

With finer frequency resolution this changes, as shown in the next figure.
 This may (as in the example) introduce additional leakage.


Matched filter - do I need to know the signal I am trying to find? I need to identify a matched filter, and so have gone through the mathematics to do this as listed here https://en.wikipedia.org/wiki/Matched_filter which derives the optimal filter $h$ to apply to a time series $x$, where that time series is assumed to comprise a desirable signal $s$ and some noise $v$. This suggests the optimal matched filter is:
\begin{equation}
    h = \frac{1}{\sqrt{s^{\mathrm{H}}\mathbf{\mathrm{R}}_v^{-1}s}} \mathbf{\mathrm{R}}_v^{-1}s
\end{equation}
and this should be applied to produce the filtered series $y$ thus:
\begin{equation}
    y\left(n\right) = \sum_{k=-\infty}^{\infty} h\left(n-k\right)x\left(k\right)
\end{equation}
however, I don't know $s$.  How, therefore, do I use this method (if at all)?

As stated by others, the signal $s[n]$ has to be known to calculate the corresponding matched filter.
An alternative is to use the power of deep learning to emulate a matched filter response as shown in this paper.

When the signal is unknown but has known statistics (PDF, covraiance etc), you match  with the estimation of the signal based on the received data, this is called the "estimator correlator" in literature. Refer to Stephen M Kay Detection theory for a more detailed explanation. The answer will be too long if I explain here.


"River" detection in text Over on the TeX stackexchange, we have been discussing how to detect "rivers" in paragraphs in this question.

In this context, rivers are bands of white space that result from accidental alignment of interword spaces in the text. Since this can be quite distracting to a reader bad rivers are considered to be a symptom of poor typography. An example of text with rivers is this one, where there are two rivers flowing diagonally.

There is interest in detecting these rivers automatically, so that they can be avoided (probably by manual editing of the text). Raphink is making some progress at the TeX level (which only knows of glyph positions and bounding boxes), but I feel confident that the best way to detect rivers is with some image processing (since glyph shapes are very important and not available to TeX). I have tried various ways to extract the rivers from the above image, but my simple idea of applying a small amount of ellipsoidal blurring doesn't seem to be good enough. I also tried some Radon Hough transform based filtering, but I didn't get anywhere with those either. The rivers are very visible to the feature-detection circuits of the human eye/retina/brain and somehow I would think this could be translated to some kind of filtering operation, but I am not able to make it work. Any ideas?

To be specific, I'm looking for some operation that will detect the 2 rivers in the above image, but not have too many other false positive detections.

EDIT: endolith asked why I am pursuing a image-processing-based approach given that in TeX we have access to the glyph positions, spacings, etc, and it might be much faster and more reliable to use an algorithm that examines the actual text. My reason for doing things the other way is that the shape of the glyphs can affect how noticeable a river is, and at the text level it is very difficult to consider this shape (which depends on the font, on ligaturing, etc). For an example of how the shape of the glyphs can be important, consider the following two examples, where the difference between them is that I have replaced a few glyphs with others of almost the same width, so that a text-based analysis would consider them equally good/bad. Note, however, that the rivers in the first example are much worse than in the second.

I have thought about this some more, and think that the following should be fairly stable. Note that I have limited myself to morphological operations, because these should be available in any standard image processing library.

(1) Open image with a nPix-by-1 mask, where nPix is about the vertical distance between letters

#% read image
img = rgb2gray('http://i.stack.imgur.com/4ShOW.png');

%# threshold and open with a rectangle
%# that is roughly letter sized
bwImg = img > 200; %# threshold of 200 is better than 128

opImg = imopen(bwImg,ones(13,1));

(2) Open image with a 1-by-mPix mask to eliminate whatever is too narrow to be a river.

opImg = imopen(opImg,ones(1,5));

(3) Remove horizontal "rivers and lakes" that are due to space between paragraphs, or indentation. For this, we remove all rows that are all true, and open with the nPix-by-1 mask that we know will not affect the rivers we have found previously.

To remove lakes, we can use an opening mask that is slightly larger than nPix-by-nPix.

At this step, we can also throw out everything that is too small to be a real river, i.e. everything that covers less area than (nPix+2)*(mPix+2)*4 (that will give us ~3 lines). The +2 is there because we know that all objects are at least nPix in height, and mPix in width, and we want to go a little above that.

%# horizontal river: just look for rows that are all true
opImg(all(opImg,2),:) = false;
%# open with line spacing (nPix)
opImg = imopen(opImg,ones(13,1));

%# remove lakes with nPix+2
opImg = opImg & ~imopen(opImg,ones(15,15));

%# remove small fry
opImg = bwareaopen(opImg,7*15*4);

(4) If we're interested in not only the length, but also the width of the river, we can combine distance transform with skeleton.

   dt = bwdist(~opImg);
   sk = bwmorph(opImg,'skel',inf);
   %# prune the skeleton a bit to remove branches
   sk = bwmorph(sk,'spur',7);

   riversWithWidth = dt.*sk;

(colors correspond to width of river (though color bar is off by a factor of 2)

Now you can get the approximate length of the rivers by counting the number of pixels in each connected component, and the average width by averaging their pixel values.

Here's the exact same analysis applied to the second, "no-river" image:

(Sorry, this post doesn't come with awesome demonstrations.)

If you wanted to work with information TeX already has (letters and positions), you can manually classify letters and letter pairs as "sloping" in one direction or another.  For example, "w" has SW and SE corner slopes, the "al" combo has a NW corner slope, "k" has a NE corner slope.  (Don't forget punctuation - a quote followed by a letter that fills the lower half of the glyph box establishes a nice slope; quote followed by q is particularly strong.)

Then, look for occurrences of corresponding slopes on opposite sides of a space - "w al" for a SW-to-NE river or "k T" for a NW-to-SE river.  When you find one on a line, see if a similar one occurs, appropriately shifted left or right, on the lines above/below; when you find a run of these, there's probably a river.

Also, obviously, just look for spaces stacked nearly vertically, for the plain vertical rivers.

You can get a little more sophisticated by measuring the "strength" of the slope: how much of the advance box is "empty" due to the slope and thus contributing to the river's width.  "w" is fairly small, as it has only a small corner of its advance box to contribute to the river, but "V" is very strong.  "b" is slightly stronger than "k"; the gentler curve gives a more visually-continuous river edge, making it stronger and visually wider.


How do you find the frequency and amplitude from a DFT? I generated a time history of 200 points using the equation
$$x(t) = A\cos(2\pi f t)+ B\sin(2\pi f t) $$

with $A = 1$, $B = 0.1$ and $f = 17.2 \mathrm{Hz}$ and a sample rate of 400 samples per second. The numerical Fourier transform of the time history has a modulus plot of:

The question is how to get the frequency, amplitude and phase from this plot? This is the typical case where the frequency falls between bins. As there are three unknowns I would imagine you could get the precise values from, say, the three largest points which I call premax, max and postmax. However, complex values are available for all the points so many methods are possible.

This question arose as an output from a previous question where it was suggested I ask afresh and that several answers were available.  I have also asked this question on Mathematica Stackexchange here where it was suggested I ask on Signal Processing.

Many answers are available. Exact rather than approximate answers are preferable. Some, like using a Prony series or curve fitting a sine wave, do not use the Fourier transform. These are also interesting but I would like to use the Fourier data.

If there is noise on the data then some sort of fitting is probably needed depending on the extent of the noise. All suggestions welcome. Thanks

Edit

A bit more in response to the answers which I still have to fully investigate

There is an exact analytic formula for the Fourier transform of a sampled sine wave. For my data it plots as

This plot could also be obtained by calculating in-between values for the original data. I also calculated the frequency at the maximum of the analytic curve. To 25 significant figures it is

17.22726210702718945079653

which is significantly different to value of 17.2 used in the calculation of the  time history.

Some of the answers suggest using parabolic interpolation. I show the results of this here

I have tried parabolic interpolation using the complex values of the Fourier transform and the absolute values. The absolute values show a better maximum than the complex values but they are still well off. I still have other suggestions to try.

Methods based on time domain analysis look more promising than using Fourier which is disappointing.

For a lower quality frequency estimate, but still more precise than just using the (center) frequency of the peak magnitude DFT result bin, try using parabolic interpolation (using 3 neighboring points) of the spectral peak.

You can't get the phase from a DFT magnitude-only plot.  For phase you will need an estimator that uses the complex DFT results.

For a higher quality estimate, including a phase estimate, try using windowed-Sinc interpolation of the complex DFT result components, perhaps followed by successive approximation to find a better estimate of the magnitude peak near the peak DFT bin.  The quality of the estimate will vary with the width of the window.

Equivalent to Sinc interpolation would be to use a much longer zero-padded DFT or FFT to get more closely spaced frequency points.

There are a few added comments regarding the above frequency estimation methods on my blog:
http://www.nicholson.com/rhn/dsp.html#1

Note that near DC or Fs/2 (closer than the width of your interpolation kernel), you may need to use a mirrored interpolation kernel.  Note that hybrid methods are possible: e.g. using Sinc interpolation to up the resolution by 4X or 10X, then using parabolic or higher-degree polynomial interpolation to estimate the peak location, then using Sinc interpolation to estimate the magnitude and phase at that estimated peak frequency.  etc. Note that sinusoidal phase estimation in noise will often be more accurate when the phase is referenced to the center of the DFT window (N/2), not to sample point 0.  Referencing phase to the center of a data vector can be done using some form of fftShift.

The question does not ask about frequency estimators that do not involve a DFT, but note that in zero noise, only 3 or 4 unaliased points are required to solve for 3 unknowns (frequency, magnitude and phase of a pure unmodulated sinusoid).  See: http://www.claysturner.com/dsp/3pointfrequency.pdf and http://www.claysturner.com/dsp/4pointfrequency.pdf

Finding the peak frequency

Your frequency resolution $\Delta f$ is directly related to the amount of samples and the aquisition frequency like this: $\Delta f = \frac{f_A}{N}$ with $N$ being the amount of samples of your signal, and $f_A$ being the sample aquisition frequency.
Thus, to increase your frequency resolution you can either decrease the aquisition frequency (Make sure to not drop below the nyquist frequency of your signal though), or increase the amount of samples.
There are 2 ways to accomplish the latter:

Increase Aquisition Time
This is pretty obvious. If you observe the signal for a longer period, you get more information about the signal, but aquiring more samples isn't always possible.
Zero-Padding
By adding zeros to your signal you can artifically elongate your aquisition time, thus generating more samples and lowering $\Delta f$. Zero-padding your signal doesn't change the information content, and is thus the most precise way to find the peak frequency.

Phase and Amplitude

Your image shows a magnitude plot of your DFT. This means that for each complex value $Y_k$ of the DFT only $\left|Y_k\right|$ is plotted. Assuming that you actually do get complex values in mathematica, and this is just the way the plotting method operates, here's what you do:
You can obtain the amplitude of a given frequency-value $Y_k$ by taking the absolute $\left|Y_k\right|$. To get the phase, you get the complex argument $\arg{Y_k}$.


Difference Between Dirac delta and unit impulse function The definition of Dirac delta function states that it gives a value of $\infty$ at t=0 and 0 elsewhere.

But, the definition of unit impulse function states that it gives a value of 1 at t=0 and 0 elsewhere.

Many textbooks state that Dirac delta and unit impulse are the same function.

So What is the difference between Dirac delta and unit impulse function?

The Dirac delta $\delta(x)$ is a continuous argument (generalized) function where $x$ is a continuous variable and can be time, space, frequency etc. Unfortunately, a simple definition like

$$\delta(x) =
\begin{cases}
\infty ~~~&,~~~ x=0 \\
0 ~~~&,~~~ x \neq 0
\end{cases}
$$

does or can not properly describe its behaviour. You need to define its properties under the integral sign, and that requires advanced mathematical stuff for a formal verification. But  a simplified Riemann definition of $\delta(x)$ is given as:

$$\int_{-\infty}^{\infty}f(x)\delta(x-a) dx = f(a)$$ for a sufficiently smooth test function $f(x)$. All properties of $\delta(x)$ follows from this Riemann (seeming) integral operation.

On the other hand, the discrete counterpart, the unit-impulse or the Kronecker delta is a sequence of integer index $m$ with very well defined behaviour as:

$$\delta[m] = \begin{cases}
1 ~~~&,~~~ m=0 \\
0 ~~~&,~~~ m \neq 0
\end{cases}
$$
These two functions serve the same purpose for discrete and continuous systems, hence share the same name impulse. Yet they are two different, distinct entities; the Dirac impulse is a very problematic concept, whereas the discrete  impulse is a mere simplicity.

Dirac Delta function(better distribution) and unit impulse are exactly same but we interpret the definition in a different way.

When we refer to this function as Dirac Delta we consider the magnitude of the distribution at every point. Hence this gives value infinite at t=0 and zero other wise.

When we refer to signal as unit impulse the values we refer are the areas of integration around that particular point with limits t- to t+. This gives value 1 to the unit impulse representation of the function.

To clarify further they are in CT domain ie continues time domain but there DT version is kroneker delta, this is defined on only integral values of n, whereas you see the unit impulse or Dirac Delta function as you may know is defined for all real values of t.

Source:Signals and Systems by  Oppenheim (Professor at MIT)


3dB and 6dB Per Octave Rolloff Digital Filters / Generating Pink and Brown Noise So I would like to have a design method for filters with 3dB and 6dB per octave roll off -- for generating pink and brown noise respectively.

I know the following 'pinking' filters exist:
Filter to add 3dB per octave?

But the poster says he doesn't remember the value for A, nor does he give any idea of how he did it in the first place.

So my main question is:

How do we pick the poles and zeros of such a filter?

-6 dB / octave is trivial.

a common reference for a pinking filter is http://www.firstpr.com.au/dsp/pink-noise/

So I have figured out the answer, thanks to a little prodding by RBJ.

Creating -6dB/octave Filters and Brown Noise:

Generate a White Noise Sequence.
Design a 1st Order Butterworth Filter
Apply 1st Order Buttworth to the White Noise Sequence
[Use the Bilinear Transform to convert the analog coefficients to digital, and apply the digital filter to the white noise sequence].

Creating -3dB/octave Filters and Pink Noise:

Generate a White Noise Sequence
Generate arbitrary interleaved poles and zeros on the positive real line within the unit circle. Such that $|p_{i}| > |z_{i}|$ and $|z_{i}| > |p_{i+1}|$. Set the first pole close to the unit circle.
Convert these roots to coefficients, and use these coefficients as the coefficients of a digital filter.
Apply said digital filter to the White Noise Sequence.

What I am still not sure of is why we need to interleave the poles and zeros, but that isn't exactly what the question asked. Can anybody elaborate?


What's the minimum decibel value? I'm trying to scale wave amplitude (0-1) to perceived loudness (still 0-1). The decibel scale seems to be unbounded. For my audio, the amplitude 1 represents a decibel $d$ - i.e. the audio is being played at volume $d$. But what would the decibel be for 0 amplitude? I know 0dB is the minimum a human can hear, but what's the actual minimum. If I can obtain this, then I can successfully convert wave amplitude to perceived loudness (while staying between 0 and 1)

By any reasonable metric, the dB value for zero amplitude is minus infinity.

As a physically realizable metric, the dB value for zero is anything less than Plank energy, or 1.96×109J, as the would be statistically unmeasurable.


Why can't a causal digital filter have an infinitely sharp transition between the passband and the stopband? In DSP book by Proakis and as well as in this pdf, it is mentioned that practical causal digital filters cannot have an infinitely sharp transition from Pass-band to Stop-band. Why is it so? Can you please provide a detailed explanation (with proof)?
Edit 1: In the book in was just mentioned that is a consequence of the Gibbs phenomenon, which results from the truncation of h(n) to achieve causality. I didn't understand how.

I don't have a concrete proof for this one. However, I can tell you this...
Consider a perfect low pass filter. The time domain representation is a sinc. And for any system to have a sharp transition band, a base signal has to be multiplied with a rectangular waveform in the frequency domain. Which implies that, the time domain signal of the same has to be convoluted with a sinc in time domain. We know that sinc function is a non causal signal.
Hence proved.

Whenever there's sharp transition in any system(filter),it'll have infinite length impulse response.For a system to be causal impulse response must be zero for negative input
h[n]=0,for n<0 ,but due to sharp transition this condition will not be satisfied,so resulting system will be non causal and non causal systems aren't realizable.


Lowest frequency within pitch detection in polyponic sound I perform some tests where I run cepstrum pitch detection algorithm over polyphonic sound. Result shows that there is almost always occurrence of lowest frequency of tones played simultaneously.
I think this is because the lowest frequency has the longest period. But I do not know if it is true or why it is so.
Can you somebody explain this or why there is lowest frequency? Or suggest some publications?

The lowest pitch in a chord is likely to have a large number of its possible overtone frequencies match the actual overtones of other notes in a harmonious chord, leading to a good cepstral magnitude, which looks for spectral periodicity (overtone series's).

However a cepstrum might also find more overtone sequence periodicities at pitches lower than the lowest pitch in a chord.


What would be the outcome of the convolution of two images? Lena

Mask

I have written a C# code which produces the following output.

Is my output correct?
Am I padding in the correct way?

You can use MATLAB to get the 2D convolution of these images. Before you do so, you'll want to convert image 2 to grayscale and convert both to double.

The result has a very large range, so I normalized the second image.

lena = im2double(imread('lena.jpg'));
color = im2double(rgb2gray(imread('color.jpg')));
out = conv2(lena, color./sum(sum(color)));

If you happen to be interested in blending images together, you might want to check out this paper.

In addition to Brian's response, this is how your "convolution kernel" looks like in "3D" (i.e, treating the colour image as a height map ):

To reproduce (Octave):

U = imread("pie.png"); %Where pie.png is a grayscale image
surf(double(U));

The reason for the blue being at the lowest level is because it has the smallest contribution to the grayscale value and conversely green is highest because it has the highest contribution.

The jagged rim is the thin border marking the circle around the hue wheel. If we now take the vertical cross section, it looks like this:

Via

plot(U(:,128));

And the horizontal cross section, looks like this:

Via

plot(U(128,:));

The spikes certainly make things interesting (because they require a lot of harmonics for their formation) but the frequency domain of such "pulses" will look a little something like this:

Via

semilogy(abs(fft(double(U(128,:))./sum(double(U(128,:))))));

And its three dimensional version looks like this:

Via

surf(fftshift(log10(1+abs(fft2(double(U))))));

(Please note, the log10(1+something) is this logarithmic transformation which previously was handled by the semilogy plot).

And if you look at it from above (just like we do with images), it looks like this:

Via

P = abs(fft2(double(U));
V = fftshift(log10(1+P));
E = V-min(min(V)); %Align the minimum to zero
E = E./max(max(E)); %Divide by maximum to align the maximum value to 1
surf(E); %And once you have the plot, rotate it to the top view, or change its view via one of the relevant functions.

Where the "varying degrees of bass" are visible (those red streaks), as we "sweep" around the complex plane at some angle $\theta$.

Hope this helps.


Invertibility of an ideal differentiator Is the system $y(t)= dx(t)/dt$  invertible or not?
If yes, please determine the inverse system for it.

It is not. It suffices to find a counter example. Let me tell you that "I am flat". Could you derive my actual value?

So, any constant function $f(t) =c$ is differentiable, and have the same derivative, $f'(t)=0$. Only from knowing that  $f'(t)=0$, you cannot recover the original (constant) function. Hence the system is not invertible.

Delta coding (with no quantization) is quite close to what you are asking. Let:

x=randi([0 10], 8,1);
delta = diff([0; x]);
x_hat = cumsum(delta);
assert(all(x-x_hat)==0);

But then you have a «silent agreement» to preprepend a zero such that you have a known starting point. And differentiation/integration is replaced by its discrete equivalents.

Knut


Echo cancellation supporting long delays and without frequency domain processing AEC algorithms mostly rely on LMS adaptive filtering, i.e. you update FIR filter coefficients then perform the filtering. Theoretically, the FIR must be as long as the maximum echo length you want to cancel. For instance to cancel delays up to 500ms on a 48kHz signal, you'll need a 24000 point FIR.
When your memory and processing power limitations make it so you can neither afford to perform 24000 MACs per processed sample, nor to use FFT-based fast convolution algorithms, is there a way around for cancelling such potentially high length echoes in a more affordable way, given delay is unknown and potentially variable?

I was wondering maybe about some other algorithm running in parallel that could assess the approximate delay, then use an adaptive length delay line + a shorter adaptive FIR filter (up to a few hundreds taps is OK)

Does this make sense? Any other neat approach to suggest?

I believe you logic is reasonable. The fact that there are companies offering algorithms for long echo tail (even above 1 second) with reasonable performance indicates that an algorithm similar to your suggested logic can work.


GNU Radio QPSK 4 Symbols I'd like to generate a simple QPSK signal in GNU radio, but so far I haven't been able to get it to perform the way I'd like and I'm not sure if this has to do with some misunderstanding I have about how QPSK works or if I just am not familiar enough with GNU Radio.

Essentially, my conception of QPSK modulation is this - 2 bits are fed into the modulator (4 possible symbols - 00, 01, 10, 11) and, based on which combination of those 2 bits was received, one of 4 possible waveforms is generated (corresponding to +/-.707 +/- .707j). Thus, when using a QPSK modulator in GNU Radio I would expect to see, in the constellation diagram, a square with 4 points at each corner, but instead I see diamonds with thousands of points, even in the absence of any introduced noise. When I print the output of the QPSK modulator to a file sink, I also see more than 4 distinct points.

Is there a way to take a random input stream of 0's and 1's and produce a modulated signal from it, like QPSK, that has only one of four values? It is important for the simulation that I'm running that there are only 4 possible modulation values before noise is added to the system. Also, I will be sending these signals between USRP's, so if there is some way to generate a QPSK signal and transmit it on a carrier frequency between USRP's that doesn't rely on the PSK Mod block (I'm not too familiar with how to implement a raised cosine filter) that would be appreciated.

It is important for the simulation that I'm running that there are only 4 possible modulation values before noise is added to the system.

So, you made a mistake when building your simulation, or your simulation already includes a frequency offset. The way you ask this, you probably did the former, but since we have zero knowledge of what exactly you've built: Can't help you with that, make sure your receiver simulation does exactly the same as what your transmitter simulation does. ("fix your program!")

I'll answer as if this was a system that's generally correctly designed:

  Thus, when using a QPSK modulator in GNU Radio I would expect to see, in the constellation diagram, a square with 4 points at each corner, but instead I see diamonds with thousands of points, even in the absence of any introduced noise.

That's not surprising! No two oscillators are exactly the same, and so your receiver and transmitter have a frequency offset.

A frequency is the linear change of phase over time. That means that your symbol phase at the receiver has a phase offset that grows constantly every time you look. Hence, the "diamond".

You need what you've been taught in your digital communications course: a frequency synchronization!

You can do that with a second order costas loop.

Hi Zman147 this is John BG jgb2012@sky.com
1. what board are you using?
2. If you haven't chosen yet and are on budget have a look at the following link

https://osmocom.org/projects/rtl-sdr/wiki/Rtl-sdr

Note that there is already a Simulink package available, so you can easily play with IQ constellations.

If you have more budget, they you should consider something like:

http://www.ni.com/en-gb/shop/select/usrp-software-defined-radio-device?storeId=11212&userLocale=en-GB&requestLanguage=en&localeGroup=en-gb&requestLangId=-1&langId=-1

hope it helps.


Higher order spline interpolation I noticed that spline interpolation with a degree higher than 3 (everything beyond cubic splines) have a very high interpolation error, hence the prediction is mostly awful. I've come across various lecture notes, slides and Youtube videos that simply indicate that cubic splines (3rd degree) are optimal and that anything beyond that is a bad idea. These sources however never mention why this is the case.

Can anyone explain to me why this is the case and maybe give me a title/link to a journal/conference paper that explains this or maybe even gives a proof.

There is no such proof because it's not always true.  It's a rule of thumb, because I guarantee that you could come up with a situation- an infinite number of situations actually- where higher order splines would do better than cubic splines.  The optimal spline order for any given situation is the exact same order as the system you are trying to model.  If the order is the same and your data points are error free (never the case, of course, except for in theoretical problems), then you should be able to model the system perfectly.

The reason they recommend not to go higher than cubic splines is because overfitting is really, really bad.  Overfitting can greatly magnify errors, while "underfitting" (choosing a spline method with an order lower than the order of the system you are modelling) introduces some low pass filtering that is either not that bad or sometimes even beneficial.

Another way to look at this is by noticing that most interesting images and data sets are "smooth" to a certain degree. Interpolating gaussian noise, for example, should work better with higher order splines.


Relationship between windowed fourier transform and wavelet transform I was reading on windowed fourier transform and wavelet transform, and i was thinking that the windowed fourier transform is a subset of wavelet transform. Is that true?

Wavelets commonly factor the frequency of the basis vectors into their temporal widths, some with more than one temporal location.  Windowed FFTs fix the temporal width of each basis vector to the length of the FFT, thus independent of any bin frequency, and with the basis location fixed to the center of the FFT.

I suppose one could consider the FFT to be similar to a wavelet transform with a bad set of wavelets (too many high frequency ones and with poor temporal/locality resolution, and too few low ones and with poor frequency resolution).

basic idea of wavelet transform is to choose basis for wavelet tranform,while for  windowed fourier transform basic idea is to choose  window  so that keeps good  frequency and time resolution,there is link about it

http://wwwmayr.informatik.tu-muenchen.de/konferenzen/Jass05/courses/2/Yakovlev/Yakovlev_paper.pdf


Time-varying "impulse response" According to many references [1,2], the time-varying "impulse response" can compute wireless channel output $y(t)$ at time $t$ using the following expression:
$$
y(t) = \int h(\tau, t) x(t - \tau) d\tau
$$

In both references, they state that this represents the response of the channel at time $t$ to an impulse applied at time $t-\tau$.

It seems reasonable to assume that there is some version of x(t) that involves a delta function that we can apply as an input that returns $h(\tau,t)$ as the output.

Trying:
$$x(t) = \delta(t)$$
$$\implies y(t) = \int h(\tau, t) \delta(t - \tau) d\tau = h(t,t) \qquad \text{nope}
$$

Trying:
$$x(t) = \delta(t-\tau')$$
$$\implies y(t) = \int h(\tau, t) \delta(t - \tau - \tau') d\tau = h(t-\tau',t) \qquad \text{nope, but closer}
$$

Is there a way to generate something resembling $h(\tau,t)$ as the output?

$$
$$

References:

[1] Proakis, Digital Communications, 5th ed, p.832

[2] Goldsmith, Wireless Communications, 1st ed, p.67

In writing $h(\tau,t)$ with $t$ is time and $\tau$ is delay, we are in the model that $\tau$ varies "differently" from $t$. Or in other words, they are different notions in spite of the fact that they are both time unit.

Similarly, in this equation
$$y(t) = \int h(\tau, t) x(t - \tau) d\tau$$

The time index of $y(t)$ varies "differently" from the time index of $h(\tau,t)$.

It is easier to write $h(\tau, t')$ and fix the measuring moment $t'=t_0$ for $h(\tau,t_0)$ and applying an input $\delta(t)$

$$y(t) = \int h(\tau, t_0) \delta(t - \tau) d\tau = h(t,t_0)$$

The index $t$ of $h(t,t_0)$ plays the role of delay $\tau$ that is exactly the same expression $h(\tau,t_0)$.

Redo the experiment at $t'=t_1, t_2, ..., $ we have the two dimensional time varying impulse reponse $h(\tau,t)$.

so, for reference (and to remind myself), the LTI convolution is:

$$\begin{align}
y(t) &= h(t) \circledast x(t) \\
     &= \int_{-\infty}^{\infty} h(t-u) \, x(u) \, du \\
     &= \int_{-\infty}^{\infty} x(t-u) \, h(u) \, du \\
 \end{align}$$

if $x(t) = \delta(t)$ then the output is $y(t)=h(t)$.

if $x(t) = \delta(t-\tau)$ then the output is $y(t)=h(t-\tau)$.

i'm gonna generalize the first integral to linear time-varying, extending $h(t)$ to two arguments.  so

$$\begin{align}
h(t)  &\leftarrow h(0,t) \\
h(t-\tau)  &\leftarrow h(0,t-\tau) = h(\tau,t) \qquad \text{(LTI)} \\
\end{align}$$

$h(\tau,t)$ is the system response to an impulse applied at time $\tau$. so, the more general convolution integral is

$$\begin{align}
y(t) &= \int_{-\infty}^{\infty} h(u,t) \, x(u) \, du \\
 \end{align}$$

so, i think your first equation is wrong.


What Is the Walsh Hadamard Transform? How Could One Use It in Image Processing? I am trying to teach myself about the WHT but there dont seem to be many good explanations of it online anywhere. I think I have figured out how to calculate the WHT, but I am really trying to understand why it is considered useful within the image recognition domain.

What is so special about it, and what properties does it bring out in a signal that would not show up on classical Fourier transforms, or other wavelet transforms? Why is it useful for object recognition as pointed out here?

NASA used to use the Hadamard transform as a basis for compressing photographs from interplanetary probes during the 1960's and early '70s.  Hadamard is a computationally simpler substitute for the Fourier transform, since it requires no multiplication or division operations (all factors are plus or minus one).  Multiply and divide operations were extremely time intensive on the small computers used on board those spacecraft, so avoiding them was beneficial both in terms of compute time and energy consumption.  But since the development of faster computers incorporating single-cycle multipliers, and perfection of newer algorithms such as the Fast Fourier Transform, as well as the development of JPEG, MPEG, and other image compression, I believe Hadamard has fallen out of use.  However, I understand it may be staging a comeback for use in quantum computing.
(NASA use is from an old article in NASA Tech Briefs; exact attribution unavailable.)

Some links:
Web Page

General description

For Gaussian distribution

Report


Problem designing a specific filter I have the next problem.

$H_{c1}(j\omega )$ is the ideal antialising filter and $H_{c2}(j\omega )$ is a real one. I'm asked to design $H(e^{j\Omega })$ so that $y[n]$ in the second diagram (the one on the right) is exactly the same as $x[n]$ in the first diagram (the one on the left).

I did something but I'm pretty sure that it isn't right. I found the equation for $H_{c2}(j\omega )$ as a function of $\omega$. Then I normalized the $\omega$-axis due to the A/D converter and divided the magnitude by $T$. If I didn't make any mistakes while doing the algebra, we get that

$$H_{c2}(e^{j\Omega })=\frac{9}{10}\cdot \frac{1}{\pi - \omega _{p}T}\cdot \Omega + \frac{9}{10}\cdot \frac{1}{\pi - \omega _{p}T}\cdot \pi + \frac{1}{10}\hspace{0.5cm} for\hspace{0.5cm} -\pi < \Omega < -\omega _{p}T$$

After that, I just thought of finding the inverse function of $H_{c2}(e^{j\Omega })$ for $-\pi < \Omega < -\omega _{p}T$... And that was it. Let $B(\Omega)=\frac{1}{\frac{9}{10}\cdot \frac{1}{\pi - \omega _{p}T}\cdot \Omega + \frac{9}{10}\cdot \frac{1}{\pi - \omega _{p}T}\cdot \pi + \frac{1}{10}}$. Then,

$$H(e^{j\Omega }) =
\left\{
 \begin{array}{l1}
  B(\Omega )  & \mbox{if } -\pi <\Omega <-\omega _{p}T \\
  1 & \mbox{if } -\omega _{p}T < \Omega<\omega _{p}T \\
B(-\Omega) & \mbox{if } \omega _{p}T < \Omega<\pi
 \end{array}
\right.$$

The problem with this is that I'm pretty sure that the filter I "designed" has no antitransform. I think that there is a smarter way of approaching this exercise but I just can't see it.

so that $y[n]$ in the second diagram (the one on the right) is exactly the same as $x[n]$ in the first diagram (the one on the left).

Considering that $H_{c1}$ is per your very own definition $\ne H_{c2}$ means that without restricting your system's input $r_c$, there's nothing you can do. You will have to implement a filter with transfer function $H_{c1}$ if you need a system with transfer function  $H_{c1}$; there's no way around that.

Now,  $H_{c1}$ being impossible with finitely long filters (as it is infinitely sharp in frequency domain), this is impossible.

I just found this exact problem in Oppenheim-Schafer's Discrete Time Signal Processing (2nd edition). For those interested, it's exercise 4.56.

I've found in Internet the solutions for those problems. The solution for 4.56 that the book gives us is

So the definition of the filter written in the original question was correct.


What is the amplitude of the impulse resonse I'm learning how to construct a discrete signal as a sum of impulse responses of the system. The final formula is then:

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

My question is, what is the amplitude of the impulse response, because i thought an impulse response was already a scaled response to the unit impulse inputted into the system? But here they are multiplying again with a scaling factor of x[n]? So is the unit impulse's amplitude just 1 then?

In case of discrete signals, the amplitude of unit impulse is 1.

It does not mean that your impulse response should have unit amplitude. The impulse response can have any real valued amplitude in it.

Note: The area under unit impulse in continuous time domain is 1 not the amplitude. The amplitude of a unit impulse approaches infinity as the area under remains constant (=1).

Impulse response means the output of an LTI system when the input is a unit impulse.

eg. $ y[n]=x[n]+2x[n-1]-x[n-2]$

is a system. where, y[n] is the output and x[n] is the input.
Then the unit impulse response of the system is calculated by substituting
$x[n]=\delta[n]$

then, the output $y[n]=\delta[n]+2\delta[n-1]-\delta[n-2]$

ie. impulse response, h[n]={1,2,-1}

Hence values in the impulse response may take any real number.

Impulse response is used to characterize an LTI system. That is if we know the impulse response of a system we can calculate the output for any input x[n] by finding convolution of impulse response (h[n]) with the input(x[n]).

Output of a system with impulse response h[n] and input x[n] is given by the equation (convolution),
$ y[n]=  \sum_{k= -\infty}^\infty x[k] h[n-k] $

Here the input x[n] does not have any connection with impulse response h[n].

Note : Unit Impulse is a standard signal with unit amplitude at 0 and 0 at all other instants.
But a unit impulse response is the output of an LTI system when the input is a unit impulse. It is the characteristics of an LTI system and the values may take any real number and is different for different system.


UWB radar to find heartbeat I have a UWB radar that Im playing with. The problem is that my algorithms at the moment can clearly, reliably and validly pick up 1-2mm movements of a pendulum i've set up. However, when I use it on my chest to find my heart, things doesnt work. The radar should be set up so that the waves penetrate skin and bone, but I am not entirely sure. The antennae are specced specifically to be able to pick up a heartbeat (or penetrate human tissues).

In terms of math, I first do despiking and a bandpass filter around 0.1 - 5 hz, followed by a principle component analysis to reduce the number of samples/sensors. Then I either do a FFT plus some stuff to find the most pronounced frequencies of the top 5 principle components, or I search for the best fit sinuswave (as the heartbeat can be described as such afaik) to the top 5 principle components.

While these methods give similar results when used on a pendulum, they break down when trying to find the heartbeat. The radar can nicely penetrate wood, and plastic and other stuff. Any ideas?

Bonus question: We had some ultrasound gel at the lab, and i tried to see if removing the "air" component would help, but then all signal dropped. Is ultrasound gel absorbing the radar pulses?

Thanks alot!

I am also working on UWB bio radar to detect the heart beats. But the problem is that breath signal is also mixed up with heart signal. I am trying to remove the breath signal (0.25 to 0.7 Hz) and to get the heart frequency signal (1  to 4 Hz). Currently I am also using an ultrasound sensor to detect heart signal.


Is there any tool to anlyze audio specturm of a song and give the average visually ? I have a lot of songs I made over the 20 years and was thinking it would be interesting to see how the overall EQ changed on the various songs I made.

Any tool out there to help me with something like this ?

Audacity will let you look at a song's spectrum/spectrogram, and it works with most audio formats.

It's also free!


How to proceess Image with curved corners and edges to give sharp Corners? I have several Images of shape H, F, I, L.. which have curved corners. How do I process them to give sharp right-angled corners? I tried approximating them with inner rectangles, but this takes a LOT of time.. How do I proceed with it?

This is the sample Image. I want to process it such that I get rid of rounded corneres.

here's how houghLinesP would look like. Not a solution yet, but maybe you can work with that:

    int main()
    {

        cv::Mat image = cv::imread("SharpenCorners.jpg");

        cv::imshow("image", image);

        cv::Mat gray;
        cv::cvtColor(image,gray,CV_BGR2GRAY);

        cv::imshow("gray", gray);

        cv::Mat output; image.copyTo(output);

        // gradient:
        cv::Mat gx, gy;
        cv::Sobel(gray, gx, CV_32F, 1, 0, 3 );
        cv::Sobel(gray, gy, CV_32F, 0, 1, 3 );
        cv::Mat gm;
        cv::magnitude(gx,gy,gm);

        //cv::Mat g_thres = gray < 100;
        cv::Mat g_thres = gm > 100;

        //int lineThres = 5;    // 1st image
        //int lineThres = 10;// 2nd image
        //int lineThres = 25;// 3rd image
        int lineThres = 50;// 4th image

        std::vector<cv::Vec4i> lines;
        cv::HoughLinesP( g_thres, lines, 1, CV_PI/(4*180.0), lineThres, 10, 10 );

        for( size_t i = 0; i < lines.size(); i++ )
        {
            cv::line( output, cv::Point(lines[i][0], lines[i][1]),
                    cv::Point(lines[i][2], lines[i][3]), cv::Scalar(155,255,155), 2, 8 );
        }

        cv::imwrite("houghpSharpen4.png", output);
        cv::imshow("g thres", g_thres);
        cv::namedWindow("output"); cv::imshow("output", output);
        cv::waitKey(-1);

        return 0;
    }

with these results for given thresholds:

1st:

2nd:

3rd:

4th:

Additional idea:

Giving your own comment, that you might try to approximate each contour with rectangles, you could start with houghLinesP results (with a pretty high threshold to get at least all "long" line elements) and try to create the rectangles based on those lines

It is easy to find the contour for each character. You can analyze the contour to find the points at the corners. Using one line to connect two adjacent points can produce right angle corners.


Frequency response creates non-causal impulse response I have frequency response data. It is output from an acoustic simulator.

I convert this into an impulse response, and get:

In that picture the blue line is the impulse response. Mostly it is hiding behind the green line (which is the impulse response if I halve the frequency resolution).

I'm sorry for an unclear picture. I will replace it in due course.

Anyway, the thing I notice of interest is that if the impulse response is regarded as a loop, the ends appear to match up.

Is there any mathematical basis for this?

It seems strange to me that the impulse response is somehow non-causal.  Because what it is modelling is a single pressure wave hitting a microphone.  Commonsense suggests to me that it should have a very sharp attack as the main wave hits followed by diminishing ripples corresponding to echoes.

So I'm quite surprised to see ripples before the main spike, let alone before t=0, as it were.

Is this some artefact of the FFT due to the fact that it is forced to sample at discrete intervals, therefore the spike cannot be completely accurate, and maybe the preceding wobbles are compensating for this.

In which case would I do better to reposition the start of the impulse response 100 samples in from the end, and have it wrap round?

I am familiar with similar phenomena. During measurements with usage of sweep-sine, when you convolve your signal with inverse filter you are obtaining your IR. Although if you consider perfect case of convolving your original sweep $x(t)$ with the inverse filter $f(t)$, then you should get perfect Kronecker delta. Unfortunately what you will get is an impulse with pre and post-ringing:

One way to minimise this effect is to use so called Kirkeby Algorithm, described here by A.Farina (my guru of acoustic measurements): Advancements in Impulse Response Measurements by Sine Sweeps. I am sure will find it quite fascinating.

One last tittle-tattle:

  Commonsense suggests to me that it should have a very sharp attack as the main wave hits followed by diminishing ripples corresponding to echoes.

You are definitely right, although in real world measurements it is possible to have preceding spike before main impulse in some cases. For example if you place your source close to the wall, then wave from speaker will excite mechanical wave in a wall and wave will propagate in the material with speed of sound much higher than in case of air. Nextly wall will emit acoustic wave that will be captured by your microphone.

impulseResponses = [impulseResponses(:,end/2+1:end), impulseResponses(:,1:end/2)]

this is the solution,
and ya remember to use a higher number of samples for fft bins,


Trignometric Fourier series representation of a continous time signal While learning Fourier series I read the definitions of representation for a continuous time signal $x(t)$ as:

$$x(t)=A_0 + 2 \sum_{k=1}^{\infty} A_k \cos(k \omega_0 t) - B_k \sin(k \omega_0 t) \tag{1}$$

where $A_k$ and $B_k$ are real.

Another definition was:
$$x(t)=a_0 + \sum_{k=1}^{\infty} a_k \cos(k \omega_0 t) + b_k \sin(k \omega_0 t) \tag{2}$$

Why two different definitions?  I mean in $(1)$ there is a coefficient of $2$ for sinusoids and sine coefficients are negated.

alright, let's review a little bit of Euler before we get to the Fourier.

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

from that you can get

$$ \cos(\theta) = \frac{e^{j \theta} + e^{-j \theta}}{2} \quad\quad\quad \sin(\theta) = \frac{e^{j \theta} - e^{-j \theta}}{2j} $$

so now let's look at Eq (1)

$$ \begin{align}
x(t) \ &= \ A_0 \ + \ 2 \sum_{k=1}^{\infty} A_k \cos(k \omega_0 t) - B_k \sin(k \omega_0 t) \\
 &= \ A_0 \ + \ 2 \sum_{k=1}^{\infty} A_k \frac{e^{j k \omega_0 t} + e^{-j k \omega_0 t}}{2} \ - \ B_k \frac{e^{j k \omega_0 t} - e^{-j k \omega_0 t}}{2j} \\
 &= \ A_0 \ + \ \sum_{k=1}^{\infty} (A_k + jB_k) e^{j k \omega_0 t} \ + \ (A_k - jB_k) e^{-j k \omega_0 t} \\
 &= \ \sum_{k=-\infty}^{\infty} c_k \ e^{j k \omega_0 t} \\
\end{align} $$

where

$$ c_k =
\begin{cases}
A_{-k} - jB_{-k}, & \quad \text{for } k < 0 \\
A_0, & \quad\quad k = 0 \\
A_k + jB_k, & \quad\quad k > 0
\end{cases} $$

and, going in the other direction,

$$ \begin{array}{lcl}
A_0 & = & c_0 \\
A_k & = & \Re\{c_k \} \quad\quad \text{for } k>0\\
B_k & = & \Im\{c_k \} \quad\quad\quad k>0
\end{array} $$

note that for real $A_k$ and real $B_k$, then

$$ c_{-k} = c_k^* = \text{"complex conjugate of } c_k \text{ "}$$

so, for Eq (1), the simplicity of having the real and imaginary parts of the $c_k$ coefficients be simply $A_k$ and $B_k$ (for $k \ge 0$) is the motivation for the convention of the leading "2" before the summation and for the minus sign.  it's just a convention.  this convention is useful because it's much easier to derive the coefficients $c_k$ than it is to derive the coefficients $A_k$ and $B_k$.

$$ c_k \ = \ \frac{\omega_0}{2 \pi}\int_{t_0 - \pi/\omega_0}^{t_0 + \pi/\omega_0} x(t) \ e^{-j k \omega_0 t} \ dt $$

where $ -\infty < t_0 < +\infty $ can be any convenient real value.

the convention for Eq. (2) comes about from looking at Fourier series first as a real analysis problem without the use of complex variables.  then this is a simple first statement:

$$x(t) \ = \ a_0 \ + \ \sum_{k=1}^{\infty} a_k \cos(k \omega_0 t) + b_k \sin(k \omega_0 t)$$

as you can see $a_k$ and $b_k$ are related to $A_k$ and $B_k$ in a very simple and straight forward manner.  you can derive that simple relationship.  but the formulae for getting $a_k$ and $b_k$ from $x(t)$ (and $\omega_0$) is less straight forward to derive and to express.  if fact, it's two equations, not one simpler equation.

The first represented is single sided (only positive frequency) that usually is introduced with reference to the double ended complex form of the Fourier series. The coefficients B and C are real because of how they relate to the generally complex coefficients from the negative and positive frequency components in the complex form.

Going back to the first form, the complex representation uses a sum from $-\infty$ to $+\infty$.  Typically a complex exponential is used in the sum, but $\cos + i\sin$ works just fine.  Say the sum is defined with COMPLEX coefficients using the letter $a$.

There is a relationship between the two sets of coefficients:

$B_k = 2\operatorname{Re}(a_k)$

$C_k = -2\operatorname{Im}(a_k)$

Hence the negative sign on $C$.

The 2nd representation is a valid representation but because the series is defined only in terms of positive frequency and real functions, the coefficients are calculated differently.

This answer is correct.  For anyone that doesn't understand what a singled ended representation of a Fourier series is, please read the following:

http://zone.ni.com/reference/en-XX/help/371361J-01/lvanlsconcepts/lvac_convert_twosided_power_spec_to_singlesided/

i.e. robert bristow-johnson

You can also consult the standard bearer on LTI systems, Athanasios Papoulis, "Circuits and Systems" pp 336.


Why do we use the HSV colour space so often in vision and image processing? I see the HSV colour space used all over the place: for tracking, human detection, etc...  I'm wondering, why?  What is it about this colour space that makes it better than using RGB?

The simple answer is that unlike RGB, HSV separates luma, or the image intensity, from chroma or the color information.  This is very useful in many applications.  For example, if you want to do histogram equalization of a color image, you probably want to do that only on the intensity component, and leave the color components alone.  Otherwise you will get very strange colors.

In computer vision you often want to separate color components from intensity for various reasons, such as robustness to lighting changes, or removing shadows.

Note, however, that HSV is one of many color spaces that separate color from intensity (See YCbCr, Lab, etc.).  HSV is often used simply because the code for converting between RGB and HSV is widely available and can also be easily implemented.  For example, the Image Processing Toolbox for MATLAB includes functions rgb2hsv and hsv2rgb.

HSV color model is more correlated to how humans see colored objects, as compared to RGB, YUV, Lab etc.

We see what color the object is (Hue), how much is it saturated (Saturation) and how much white light is falling on it (Intesity).


Use of pulse shaping in digital modulation Sorry, I am a little bit confused here cause I thought when we say QPSK, the signal looks like the following.

I thought the bits to send/receive depends on the phase shift of the carrier, or is the above picture wrong or will never happen in real world?
So the discussion from here is to use pulses to send 1's and 0's? What does the signal look like in air actually? Pulse Shaped or like my attached picture? I am confused, so does the phase shifted carrier signal indicate a symbol or does the amplitude of a pulse shaped signal carries the bits?

QPSK in its raw form would appear just as the OP has shown, where the information is encoded into four symbols given by the four phase states alone. The problem is this has a very wide transmit spectrum given the pulse shapes are rectangular and the Fourier Transform of a rectangular pulse is a Sinc function (whose peaks roll-off in frequency only at the rate of 1/f). Wireless spectrum is an expensive commodity, so we go through great lengths in terrestrial wireless communications to limit spectral occupancy and in the process improve spectral efficiency. This is the only reason for pulse shaping, but it is significant. By more slowly transitioning between symbols rather than the abrupt transitions shown, we greatly limit the amount of spectrum needed (the spectrum will be given by the Fourier Transform of the pulse shape). Effective pulse shaping does this while controlling the inter-symbol interference that is added in the process.
The graphic below demonstrates this for a 16QAM waveform where the real portion of the time domain waveform is shown in the upper part of the graphic with and without pulse shaping. The lower part of the graphic shows comparative spectrum where we see the significant reduction in spectral occupancy due to the pulse shaping provided.

I see this graphic is the same one as what the OP linked, so to further explain: This is how the envelope of the signal would look like in the air (we don't see the actual carrier frequency in this graphic). This graphic was for QAM so does have information in the amplitude as well as phase, but for QPSK case the amplitude would carry no information. But even with the QAM case, the pulse shaping modifies the amplitude ONLY to reduce the spectrum, and modifies it in such a way to control the trajectory between our symbol samples of interest, while still passing through those samples defined by the modulation exactly at the proper symbol timing locations.
Below is an eye diagram showing the real and complex waveform for the QPSK case (The eye-diagram similar to the time domain waveform shown above, but repeating synchronous to symbol boundaries)

Show this shows how the over the air waveform would look (in actuality we typically do half the pulse shaping in the transmitter with a root-raised cosine filter while this plot shows the result after the second root-raised cosine filter in the receiver), but in general we are seeing the amplitude and phase of the time domain waveform for all samples, and the job of the receiver is to determine the sampling locations where the red dots are shown (the center of each symbol), and we then determine the magnitude and phase of those to demodulate (or demap) the QPSK data.
To see what is transmitted "over the air", the sketch below begins to demonstrate how the envelope of the carrier would be modulated by the pulse shaping, but doesn't provide the whole picture as the phase would also be transitioning from symbol to symbol so would not align with what would occur if we continued to transmit one symbol repeatedly. I added more details further below showing the exact waveforms expected. These are the real waveforms that are transmitted at the antenna rather than the baseband equivalent waveforms showing the complex envelope typically with two plots as I (real) and Q (imaginary). The sinusoidal carrier could be any frequency and the envelope would not change, including frequency = 0 which would then represent the eye diagram and I/Q plots that I have shown above (and must then be a complex waveform to represent).

A further note is the OP used symbol mapping that is not optimally gray coded (the 180 degree transition depicted in the diagram going from symbol 00 to 10 should be reserved for transitions where both bits change such that only one bit changes for all symbols that are close to each other).
Also note how this is similar to the operation of windowing prior to computing the FFT to reduce spectral leakage.

To avoid any confusion based on subsequent questions in the comments, here is the actual results on what the OP's QPSK waveform would look like at "RF" using properly gray-coded symbols using the following mapping of symbol to constellation:
0: -1-1j
1: 1-1j
2: -1+1j
3: 1+1j
We will start with the final waveform in the receiver as the desired QPSK constellation and work back to what is transmitted over the air by adding in the carrier offset and seeing the result:
First the I and Q waveform at baseband for a QPSK waveform pulse-shaped with a raised cosine filter with alpha = 0.3 would appear as the following in the receiver once all the carrier offsets are removed:

The constellation of the above is formed by plotting I vs Q on the complex plane. These trajectories represent the phase and amplitude of our RF carrier at every moment in time, and here I show the one instance that occurs at every integer symbol number in red, which is the only point over the duration of each symbol where the phase and magnitudes of the carrier will actually represent what we expect for QPSK (every other sample in between will be in transition at other phases and amplitudes!).

With an actual carrier frequency added, this same sequence would become the waveform shown in the plot below, where the amplitude of the envelope will match the amplitude of the constellation in the plot above:

Zooming in on the first 10 samples, this is what the OP would see for the sequence given (this is for the sequence 0,2,1,3,0,2,1,3,1,1,2,2.... specifically). The symbol sample points are aligned in this plot with the integer symbol locations on the horizontal axis.  This was done with a carrier frequency that is 10 times the symbol rate, but if we increased the carrier frequency further this plot showing the envelope would just get solid but appear identically the same otherwise.

And the same thing as we reduce the carrier to 4 times the symbol rate, showing that we will see the same basic envelope shape regardless of which carrier frequency is used (and as noted if the carrier gets too small compared to the symbol rate we need to use a complex signal to properly represent the waveform, exactly as we do in transceiver implementations):

Note that these are all done with a raised cosine filtered pulse-shape, but what is typically transmitted over the air is root-raised cosine filtered.  To be complete, below is the plot showing the same waveform with a root-raised cosine pulse shape:

With close observation of the above plot against the unmodulated carrier aligned with the first symbol we see how the modulated carrier is continuously changing in phase during the course of each symbol and, as also made clear by the constellation diagram of the modulated waveform, is only at the correct amplitude and phase according to the QPSK constellation only over very short durations, and in this specific case after being passed through a second RRC filter. (If the pulse-shape is altogether eliminated then the waveform will stay at the correct amplitude and phase for the entire symbol, and as the pulse shape roll-off factor is increased, this will approach that).

The above plots and descriptions used Raised Cosine (and Root-Raised Cosine) pulse shaping but an alternate harris-Moerder filter is described in this paper http://eon.sdsu.edu/~seshagir/SDR05.pdf. Using this as the pulse shaping filter results in an order of magnitude lower EVM for the same number of taps. This isn't commonly used even today (15 years later) to my knowledge since the root-raised cosine filter is baked into our specifications but certainly is very compelling to be aware of for consideration in future radios (I added a question to see if there is further experience with its use: Use of the harris-Moerder Nyquist Pulse Shaping Filter which may have additional information on this.)

The QPSK signal never looks like your attached picture. Neither in the transceivers, nor "in the air".
First, I cannot remember an application where the QPSK symbol rate is on the order of the radio frequency. The 802.11a's symbol duration is 4 usec, with 5GHz carrier it gives 20,000 wave periods. The bandwidth of conventional narrowband and carrier wave transmissions is usually a tiny fraction of the carrier bandwidth. Otherwise, this communication is called Ultra Wide Band, a rapidly advanced technology.
This is not to say that in simulation we cannot modulate a carrier with the signal at rates that high. It may have no practical use, bit still be an instructive experience, to perform these simulations "by hand" and visually inspect the results. The two pictures below show the results of this simulation.

The first picture follows your attached picture in that there are two carrier wave periods per symbol.

The top graph is the (real part of) product of complex symbol amplitude and the carrier. Notice that the phases in QPSK are 45°, 135°, 225°, and 315°, in contrast to 0°, 90°, 180°, and 270° of those shown in your picture. I will not elaborate on the subject, only add that there are no pure carrier wave present in the transmission that one would be able to filter out from the antenna signal. The "carrier waves" at the centroid lines of the signal graphs can be considered "ghost carriers".
The middle graph is the symbol-carrier product filtered out with the sinc filter (brick filter in the frequency domain).
The bottom graph is the symbol-carrier product filtered out with the pure raised cosine filter, the roll-off factor = 1.

The second picture shows the waveforms with eight carrier wave periods per symbol, the intermediate case between a narrowband transmission and UWB.
The answers would be more more helpful and relevant, if you first try to run simulations for yourself, and ask the SE community for help in writing the simulation code, if you need, or for interpreting simulation results. Instead, with the question like this, the posters often guess what you really need to know.


Comfort Noise Insertion in VAD algorithms I am in a intention of Porting Voice Activity Detection algorithm in real time hardware,so i was successful in classifying speech and silence/noise but now the problem is replacement of silence frames with comfort noise(preferebly any stationary noise)

so can any one suggest any opensource CNG algorithms or how a comfort noise can be of any type(pink,white etc)

to do this usually we use a rand() functions but where as these type of functions are not available with real time DSP so how to port such type of Algorithms?

You may want to try a two-pronged approach.  First, have a standard pre-packaged CN to inject (like jan and Jason R suggest) at the start of the conversation/exchange.

As the conversation progresses, gather low-volume data, perhaps as part of your noise reduction algorithm I assume you already have, by sampling the lowest volume signal coming in between VAD end and start points. Compare your low-volume samples to your standard CN, and replace (or convolve to simply dilute your synthetic CN) if it matches some threshold or characteristic you established. All you need is some rand and normalization functions, and perhaps FFT, all of which are readily available in C (and assembly).

Here to my question i have gone for a simple solution as suggested by some people above and also with my idea i stayed simple with the approach

on the decoder side i have taken two types comfort noise with different noise strengths,while streaming the decoded data if the user finds distrubing with one noise strength i will switch to the other noise

i have got the noise files of 1 second up on which repetition will also doesn't bother a human ear


Can a wavelet transform give time dependent phase of sinusoids in signal I have a signal which contains sinusoidal components that oscillate at different frequencies. I think the phase of the sinusoids is changing with time. I could do a Fourier transform on small chunks of the signal and the phase of the sinusoids from that, but is there a better way to do this with a wavelet transform? I'm new to the concept of wavelets. From what I understand, it's good at doing what an FFT does but with time resolution. Could you get the phase of a sinusoid as a function of time from a wavelet transform ?

It is indeed possible (up to a period), with complex "continuous" wavelet transforms, provide your mother wavelet ooscillate sufficiently. You can look at the scalogram modulus and its phase at the appropriate scales, as illustrated from Matlab: Wavelet coherence:

Alternatively, you can convert your signal into an analytic one (using the Hilbert transform) and perform a real wavelet analysis.

Finally, you might dig into other linear (short-term Fourier transform) or bilinear time-frequency (Wgner-Ville, Choi-Williams) or time-scale tools.

I would look into the short-time fourier transform. The wavelet transform is closely related to the STFT, more so than the FFT. I am not sure about getting phase information, but if you want to understand the wavelet transform more, start with the STFT.


Computing impulse response for system with feedback terms I'm trying to compute the impulse response for the following system:

$$
y(n) = \frac{1}{2} \left( x(n) + x(n-1) \right) + \frac{1}{2} y(n-1)
$$

and am told to assume that $y(-1)=0$.

My solution: Letting $x(n) = \delta(n)$, we have

$$
h(n) = \frac{1}{2}\left(\delta(n) + \delta(n-1)\right) + \frac{1}{2}h(n-1).
$$

Plugging in a few values for $n$ (and assuming that $h(-1)=0$), I've obtained:

$$
h(0) = \frac{1}{2},
$$

$$
h(1) = \frac{1}{2} + \frac{1}{2} \cdot \frac{1}{2} = \frac{3}{4}
$$

and

$$
h(2) = \frac{1}{2} + \frac{1}{2} \cdot \frac{3}{4} = \frac{3}{8}.
$$

From this, I've observed that

$$
h(n) = \begin{cases}
0, & n<0 \\
\frac{1}{2}, & n=0 \\
\frac{1}{2} + \left(\frac{1}{2}\right)^n, & n > 0
\end{cases}
$$

Does this seem correct? I feel like I made a mistake somewhere and that the expression for $h(n)$ should be in a simpler form.

Note: I computed the impulse response by taking the Fourier Transform and then plugging in values for $n$ and got something different than this. However, I believe that this discrepancy is due to the initial condition $y(-1)=0$.

Your values for $h(0)$, $h(1)$ and $h(2)$ are correct. However, the general formula isn't.
It should be

$$h(n) = \begin{cases}
0, & n<0 \\
\frac{1}{2}, & n=0 \\
\frac{3}{2}\left(\frac{1}{2}\right)^n, & n > 0
\end{cases}$$

A simple way to solve such problems is to use the $\mathcal{Z}$-transform:

$$Y(z)=\frac{1}{2}X(z)(1+z^{-1})+\frac{1}{2}Y(z)z^{-1}$$

which gives

$$H(z)=\frac{Y(z)}{X(z)}=\frac{1}{2}\frac{1+z^{-1}}{1-\frac{1}{2}z^{-1}}=
\frac{1}{2}\frac{1}{1-\frac{1}{2}z^{-1}}+\frac{1}{2}\frac{z^{-1}}{1-\frac{1}{2}z^{-1}}$$

With the transform pair

$$\frac{1}{1-az^{-1}}\Longleftrightarrow a^nu(n),\quad |a|<1$$

and with the property

$$x(n-1)\Longleftrightarrow z^{-1}X(z)$$

you arrive at

$$h(n)=\frac{1}{2}\left(\frac{1}{2}\right)^nu(n)+\frac{1}{2}\left(\frac{1}{2}\right)^{n-1}u(n-1)$$

which corresponds to above result.

The line $h[2] = \dots$ should be $h[2] = \frac{1}{2}h[1]$

You can use the z-transform to derive the inpulse response:
$x[n] = \delta[n] \Rightarrow X(z) = 1$
$Y(z) = .5 (1+z^{-1}) + .5 Y(z) z^{-1} \Rightarrow Y(z) = .5 \frac{1+z^{-1}}{1-.5z^{-1}} = .5 \left( \frac{1}{1-.5z^{-1}} + z^{-1}\frac{1}{1-.5z^{-1}} \right)$
$\Rightarrow y[n] = .5 \left( .5^{n} u[n] + .5^{n-1} u[n-1] \right)$, Where $u[n]$ is the heaviside step function.

See wikipedia for transform pairs.


Can software distinguish between different musical instruments If possible, I would like to write a little script that can identify different musical instruments from just a recording.  Apart from variations of amplitude and frequency over time (which can be shown by time-dependant and spectrogram, respectively), what other variables are there?

Furthermore, how would one do the classification?  Fingerprinting is an option, but that would only work if a single piece at a time is played.

What you want to do is timbre recognition.

I think this is a bigger task than "a little script", but could be done with some determination.  I've had the same idea over the years, and this is how I'd go about it:

You'd be working with individual tracks of reference instruments, and would build an 'overtone array' of the relative power of the overtones as compared to the fundamental frequency.  You would also need to write a script to detect the fundamental frequency, of course.  Other features would also be of interest like attack and sustain characteristics.

Once you had that information for all of your reference instruments, you would take the same measurements of a new recording and compare it to your reference overtone arrays and other characteristics, and output whichever one it is most similar to.


Incoherence: Compressed Sensing (CS) vs Matrix Completion (MC) I am seeking a clarification of the concept of Incoherence within the MC framework.

Specifically,

1) the literature mentions the application of a "strong incoherence" given a set of assumptions. Intuitively, how is the strong incoherence different than "incoherence".

2) Does the Restricted Isometry Property have anything to do with low rank Matrix Completion of natural images, or does it only pertain to the reconstruction of signals in CS (i.e. y = Ax).

Background

There is an incoherence property $\mu$ used in CS to measure the correlation between the sensing and representation bases; whereby low correlation enables signal reconstruction of sparse signals with fewer samples and high correlation makes it nearly impossible regardless of the sparsity of the signal.

For MC, there is also an incoherence property $\mu$: it measures the correlation between the singular vectors of the matrix (U and V from SVD) and the standard basis vectors (i.e. $\ e_i$ = vector of zeros except for a 1 at index i). As with CS, low rank (i.e. singular value-sparse) matrices exhibiting a low correlation are more likely to be reconstructed than those with high correlation.

Both of these incoherence properties are used to stipulate the minimum number of samples required to reconstruct a signal (1D, 2D, respectively). However, I cannot see where the matrix incoherence property addresses the nature of the sampling; I.e. it only tells me which matrices are more likely to be reconstructed. I believe this is because the "nature of the sampling" that is assumed (by CS community doing MC) for low rank MC problems is strictly "uniform random sampling".

Now, restricting our attention only to those low rank matrices with low $\mu$, is there an additional MC incoherence property that characterizes the nature of the sampling to say whether or not a matrix can be reconstructed? (I'm looking for a link back to CS where they compare the relationships between different representation/sensing bases to see if they work well together to reconstruct the signal.)

Thank you.

Partial answer:

2) Matrix completion does "not rely on any notions of restricted isometry." Candes and Recht, 2008.


Does $ Ev(\cos(4\pi t)u(t)) $ signal has period? According to Oppenheim's signal and systems solution book and other solutions like this

$$
Ev(\cos(4\pi t)u(t)) = 1/2 \times ( \cos(4\pi t)u(t) + \cos(4\pi t)u(-t)) = 1/2(\cos(4\pi t))
$$

$$
\textrm{for }-\infty < t <  \infty
$$

is periodic and its period $1/2$.

But I think it is not equal $ 1/2(\cos 4\pi t)  $ at $ t = 0 $ so it is not periodic because at $ t = 0 $ it is equal to $ 1 $ different from other multiples of period.

I have also curious about the value of $ u(t)+u(-t) $. Is it $0$, $1$, or $2$ ?

It is a matter of convention. The answer depends on the definition of the step function $u(t)$. In signal processing, it is common to define $u(0)=1/2$. Under this definition, we have that $Ev(\cos(4\pi t)u(t)) = 1/2(0.5+0.5)=1/2$, so the signal is indeed periodic.

If $u(0)$ is defined as either 0 or 1, then you'd be correct and the signal wouldn't be periodic.

x(t)=1/2(cos 4 pi t u(t)+cos 4pi(-t) u(-t)) ... even function
cos(-x) =cosx but u(-t) =o as unit step function exists when t>=0
so u(-t) =0 u(t) =1 for all even values
x(t)= 1/2 cos 4 pi t ... u(t)=1 for all even values
t= 2 pi / omega  omega = 4
t= 2pi/4 =1/2 sec [ periodic function]

if cos is replaced by sin then sin function gets 0 for even values t=0 u(t) =1 but sin 0 =0 .. x(t) =0 non periodic


Convolution by differentiation property of convolution This is the question

when i am trying to do it by the deravative property of convolution
https://en.wikipedia.org/wiki/Convolution#Differentiation like this

Actual answer should come like this

  $$y(t)=
\begin{cases}0, & \text{for }t < 0 \\
7-3t, & \text{for } 0<t<1 \\
5-t, & \text{for } 1<t<5 \\
7, & \text{for } t>5
\end{cases}$$

i am getting answer -6 at t=4 bu the actual answer should come 1

  P.S:-
  Please dont give me alternative method to solve this i know it already instead tell me whats wrong i am doing in this method
  also i know i am violating the policy by putting the image but i will edit it lator in text form .

First of all let me state that, the property of differentiation under (convolution) integral operator may not be holding true for arbitrary signals but only for some subclass, such as absolutely integrable ones. That's stated in the wikipedia link OP provides about convolution. All I can make sense from that statement is that interchange of differentiation and integration is only possible if the integral is convergent. That's formal mathematics and I'm not very sure whether that applies to generalized functions, such as $u(t)$ and $\delta(t)$ as well. So I'll leave the discussion to mathematicians for a moment and proceed with the solution, treating the functions as if they're well behaving.

So far you've correctly obtained the signal $y(t)'$ from the differentiation of convolution property as: $$y(t)' = h(t)' \star x(t)  = -3 u(t) + 2u(t-1) + u(t-5)$$

Now it comes to proper integration of $y(t)'$ to obtain $y(t)$ which was the most misleading step; $$ \int_{-\infty}^{t} y(\tau)' d\tau = \int_{-\infty}^{t} -3 u(t) + 2u(t-1) + u(t-5) d\tau$$

the left side evaluates as;

$$ \int_{-\infty}^{t} y(\tau)' d\tau = \int_{-\infty}^{t} dy d\tau = y(\tau)+C |_{-\infty}^{t} = (y(t)+C)-(y(-\infty)+C) = y(t)-y(-\infty) $$

Now here I shall make a few simple comments. As below we argue with @msm whether this is a definite integral or an indefinite one? He seems to claim it as an idefinite one with constant of integration at the output. (I don't actually know it because he doesn't put his equation there inspite of all my attempts, also he seems(?) to delete that comment as well? Rather he says it's simple to see. I don't know. It's always better to make it clear when requested to avoid any confusion. Of course this is a definite integral with clear limits on it, hence at the output you won't have a constant of integration whatsoever. The integrand is antiderivatived, and evaluated at the limits and subtracted so that the constant of integration is cancelled. What remains there as $y(-\infty)$ the integrand evaluated at the limits. For causal systems and signals applied, $y(-\infty)$ is by defintion zero, hence it's ignored, but when you have a noncausal system (like in the example system with $h(t) = u(-t-1)$ you may not ignore it and hence the term $y(-\infty)$ remains.

The right side of the equation is integrated as follows:
$$  = \int_{-\infty}^{t} -3 u(\tau) + 2u(\tau-1) + u(\tau-5) d\tau  = \int_{0}^{t} -3 u(\tau) + 2u(\tau-1) + u(\tau-5) d\tau $$

$$  = (-3 r(t) + 2r(t-1) + r(t-5) + C)-(-3 r(0) + 2r(0-1) + r(0-5) + C)$$
$$  = -3 r(t) + 2r(t-1) + r(t-5) $$

Where we have indeed antiderivated each integrand and subtracted their evaluation at $\tau=t$ and $\tau=0$, the latter one all providing $0$

Hence we have;
$$ y(t) = -3 r(t) + 2r(t-1) + r(t-5) + y(-\infty) $$

At this point it's worth noting that, the output $y(t)$ depends on itself due to the $y(-\infty)$ term. Where you can fortunately (easily) evaluate it from the convolution integral as $$y(t) = \int_{-\infty}^{\infty} h(t-\tau)x(\tau)d\tau $$ as $t \rightarrow -\infty$ resulting in;
$$ y(t) = \int_{-\infty}^{\infty} h(-\infty-\tau)x(\tau)d\tau = y(t) = \int_{-\infty}^{\infty} 1 x(\tau)d\tau = 7 $$

Hence the solution becomes, $$y(t) = 7 -3 r(t) + 2r(t-1) + r(t-5) $$
and finally $$y(4) = 7 - 3 r(4) + 2 r(3) + r(-1) = 7 - 3 \times 4 + 2 \times 3 + 1 \times 0 = 1$$

Below you can see the verification by direct evaluation of the colvolution integral.

(EDIT: see the first paragraph to get the point)
As warned by the OP, that the differentiation property under the convolution is not valid for the signal $h(t) = u(-t-1)$, so we then resort to direct evaluation to check the answer.

$$y(t) = x(t) \star h(t) = x(t) \star u(-t-1)$$

The graphical method provides the easiest evaluation of this integral for particular output value of $t=4$ as:

$$y(4) = \int_{\tau= 5}^{\tau=6} x(\tau)h(4-\tau) d\tau = \int_{\tau= 5}^{\tau=6} 1 \cdot 1 d\tau = 1 $$

So the requested value is: $y(4) = 1$ , hence you are right the answer is $y(4) = 1$ and (EDIT: see the first paragraph to get the point) [it seems the differentiation property is not applicable for arbitrary signals but those which are absolutely integrable.]

An example to see the diferentiation property is; $x(t) = u(t) - u(t-1)$ and $h(t) = u(t) - u(t-1)$  where we have $y(t)' = u(t) - 2u(t-1) + u(t-2)$ and the differentiation property yields: $$y(t)' = h(t)' \star x(t) = (\delta(t) - \delta(t-1)) \star x(t) = u(t) - 2u(t-2) + u(t-2)$$

Sorry for the confusion...


Is this system an FIR filter or an IIR filter? Consider this system:

 y[n]=y[n-1] + x[n]

It can be verified easily that this system has a finite impulse response (by putting x[n]= delta[n]  i.e impulse, and evaluating impulse response h[n]=y[n]). Yet for a unit step input U[n], this system is an unstable system.  On wikipedia it is mentioned that FIR filters are always stable

Now my question is , is this an FIR filter or an IIR filter?

It is an IIR filter.  I think you misunderstand what is meant by "finite impulse response".  It does not mean finite amplitude, it means finite time- i.e. that the impulse response eventually becomes zeros.

The impulse response of the filter is a step function, which is infinite in time.

x[-1] = 0
x[ 0] = 1
x[ 1] = 0 ...

y[-1] = 0
y[ 0] = 1
y[ 1] = 1
y[ 2] = 1 ...

It is obviously IIR because the output sequence y[...] occurs on both the left-hand and right-hand sides.


Modeling channel response as complex vs splitting I and Q If I have a communications channel where a digital QAM signal is transmitted, I can either model the channel's impulse response as a complex $h(t)$, or as 2 separate impulse responses $p(t)$ and $q(t)$ for the I and Q channels, respectively. Let $s(t)$ be the transmitted signal, with $s_I(t)$ and $s_Q(t)$ denoting its I and Q components.

For the complex case, the filtered signal would be $$[s_I*\operatorname{Re}(h)-s_Q*\operatorname{Im}(h)]+j[s_I*\operatorname{Im}(h)+s_Q*\operatorname{Re}(h)]$$

For the 2 separate case, it would be $$[s_I*p]+j[s_Q*q]$$

Clearly, they look rather different. Is either one a "better" model than the other? Better in the sense of, for example, using an adaptive filter to equalize the channel. We can train for a complex response with the cross coupling of the former, or run the channels separately. If not, can we talk about some sort of mathematical equivalence between them? Is there a mathematical mapping?

After giving it more thought, I see my original confusion. We need to be careful about what we mean by the channel's impulse response. Normally for an equalizer for this application, we are not modeling the channel alone, but rather also including the frequency mixing up and down, before and after, transmission through the channel. When we include the mixers in this way, we get the cross coupling as I mentioned for the first case. So if $h(t)$ is the response of the channel alone, then $g(t)=e^{-j\omega_ct}h(t)$ is the response of the entire system from baseband to baseband. The second case with independent responses would be for applications where the channels are not mixed together at any stage and are transmitted and received independently.


Is DSP a branch of Computational Science? The way I understand it, computational science covers a wide range of fields, all of which use computers and mathematical models to solve a problem involving one of the sciences (physical, life etc.).

But is DSP, which does use maths to manipulate a signal, a form of computational science?

Nope. Its not . Although , computational sciences covers very wide range of subjects but to be able to implement DSP , one must have a thorough understanding of signal fundamentals . That is not a part of computational sciences .

DSP is a branch of electrical engineering, not science.


DFT coefficients meaning? What "are" they? What's a sensible way to interpret the coefficients (and what isn't)? To pose specifics:

DFT coefficients describe the frequencies present in a signal
They describe the sinusoidal frequencies of the source

Are either of the above true?

At it's most fundamental, the DFT is about fitting a set of basis functions to a given set of sampled data.  The basis functions are all sinusoidal functions, expressed as the complex exponential with a purely imaginary exponent.  Using the most common scaling convention each basis function, without its scaling coefficient, is:
$$ g_k[n] \triangleq \tfrac1N e^{j \omega_k n} $$
For the DFT:
$$ \omega_k = \frac{2 \pi k}{N} $$
and the basis functions add up as
$$\begin{align}
 x[n] &= \sum\limits_{k=0}^{N-1} X[k] \ g_k[n]\\
\\
  &= \sum\limits_{k=0}^{N-1} X[k] \ \tfrac1N e^{j \omega_k n} \\
\\
  &=  \tfrac1N \sum\limits_{k=0}^{N-1} X[k] \ e^{j \frac{2 \pi k}{N} n} \\
\end{align}$$
It's real easy to solve for the coefficients:
$$\begin{align}
 X[k] &= \sum\limits_{n=0}^{N-1} x[n] \ e^{ -j \frac{2 \pi k}{N} n} \\
\end{align}$$

Note: Originally I failed to mention the "complex" clarification to follow (and at bottom), hence the downvotes (which now are undue). The coefficients give amplitude and phase of a complex sinusoid. So instead of just summing real sinusoids, we're also summing imaginary (to a nonzero).
This terminology is also used ubiquitously in Wavelet Tour, with one line reading near-identical to mine (pg 26):

(1) - no. (2) - no. But not entirely no.
The truest meaning stems from the definition of forward and inverse transforms, but it won't directly answer all pertinent questions; latter is a matter of what said definitions imply. The formulation I find most intuitive stems from the inverse:
DFT coefficients describe the sinusoids that, when added, yield the original signal. -- In detail:
DFT coefficients, $X_k$, give amplitudes and phases of sinusoids at integer frequencies $k$, from $0$ to $N-1$, that sum to the original signal $x[n]$, comprised of $N$ points.
That's all. This is the one description that must hold no matter the context. (1) and (2) "turn out true" GIVEN certain criteria are met, but the coefficients themselves don't inherently warrant any of them.

Case 1: $f=1, 4,$ and $8$ adjacent, each lasting for 100 samples.

There are only three frequencies, so if the coefficients indeed describe "what frequencies are in a signal", what's with all these other frequencies?

"There are three peaks, one for each frequency; other values are due to sharp jumps between frequencies, insufficient samples for proper representation, and other imperfections."

Nonsense. There's nothing 'imperfect' about the coefficients; they simply don't directly answer for "which frequencies are in the signal". Instead, they give the amplitudes, phases, and frequencies of sinusoids that, when added, equal the original signal.
This is sometimes formulated as "Fourier Transform tells us which frequencies are in the signal, not when they occur". True, but it misses the point.

Case 2: $f=1$, 100 samples, with a 1 sample "dip".

What's with all the higher frequencies, even if they are relatively small? Surely this is noise, and not part of the 'actual' signal, and we should just pretend they're zero - right?
Wrong. These can very well be 100 different sources overlapping exactly as described by coefficient to produce this exact dip. It could also be noise. We don't know. Of course we can make rational assumptions about what's plausible, but this requires knowledge about the source system; how many sinusoidal sources are there, if any? What noise can we expect? The example above is trivial; if we know nothing else, it's almost certainly noise - but a real-world signal is far more intricate.
This addresses both (2) and (1). Not only "we don't know", but we can't know anything exactly the source system just by observing it (its "output"). It could be a magnetic pole spun at 1Hz, with an electrical shock inserted in the circuit at an instant. Or it could be 100 antennas synched just right. The end-result is identical, and it's the only thing the DFT operates on.

Case 3: $f=1$.

Now we know the source is a sinusoid of $f=1$, right? It has to be!
Nope. The source could very well be a bunch of triangle waves, square waves, a happy coincidence, or infinite other possibilities. But what if we know it's periodic? Even then: sinusoids are only one, not only, example of periodic orthogonal basis functions.
PS, the slight non-zero values in the coefficients in this case are due to 100 samples imperfectly representing an $f=1$ signal. It's also true in Case 1, but is far besides the point.

So when do the coefficients exactly represent the "frequencies in the signal", if ever?
Only one possibility; if the signal is:

Comprised of sinusoids ...
of the exact frequencies ...
at the exact amplitudes ...
at the exact phases that are given by the coefficients, ...
each lasting over the entirety of the transformed signal, ...
and nothing else.

In other words, if it does accurately describe (1), it only "happens to be the case" (i.e. coincidence). The more the source deviates from any of the six above, the less the coefficients meet (1) or (2).
If the DFT can't tell us anything about the source, why use it? -- Because most of the time we have some knowledge of the source, allowing us to assume some traits and deduce others (e.g. electric motors). Other times we don't care about the nature of the source, but only the signal and how we can manipulate it (e.g. recorded audio, power transmission).

What "is" the DFT? -- A transform. A mapping. An algorithm. Something that takes numbers, and returns other numbers, according to some rule. What it's not is a descriptor of the "true nature" of the signal.
Does DFT assume the input's periodicity? -- Nothing about the coefficients or the reconstructed signal change regardless what the signal is outside the portion that was transformed. The DFT describes only that portion, nothing else. The extension of the inverse, however, is periodic, and this becomes relevant in certain operations (e.g. below) - but again, no info on the original signal.
But DFT derivation from the Fourier Transform is only valid if input is periodic. -- The DFT neither requires nor assumes FT's existence; it is standalone. "From FT" is only one derivation, not sole. Where periodicity is relevant is when FT's properties are assumed to also apply to DFT.
In circular convolution, multiplying FT's coefficients exactly corresponds to time-domain convolution, but not DFT's; for DFT one must pad, which changes the coefficients. So in this sense, the interpretation of DFT coefficients changes with the periodicity assumption.

What about the meaning from the "forward transform" perspective? This is more about how we get the coefficients in the first place, and why they take on the values they take, which is its own topic. For start, it's better posed as a question:
Why does the forward transform yield the coefficients that describe the frequencies, amplitudes, and phases of sinusoids that sum to the (transformed) signal?
I defer to this excellent video and this text for exploring the answer.

NOTE: Answer uses real sinusoids to illustrate concepts for simplicity, which holds for real-valued inputs. However, this isn't the whole story, and may mislead, but we can apply the same concepts to the complete picture. The more complete definition adds one word:
DFT coefficients, $X_k$, give amplitudes and phases of complex sinusoids at integer frequencies $k$, from $0$ to $N-1$, that sum to the original signal $x[n]$, comprised of $N$ points.
For a real-valued input, due to symmetry of $X$'s imaginary components, the imaginary components of (coefficient-multiplied) complex sinusoids sum to zero. However, the imaginary component of the complex sinusoid (basis function) does not contribute zero; it contributes to the real component:
$$
(A + jB) e^{j\theta} = (A + jB) (C + jD) = (AC - BD) + j(AD + BC), \\
C = \cos{(\theta)},\ D = \sin({\theta})
$$
Thus, the real sinusoids used throughout the answer actually owe part to the imaginary component of the complex basis function, so it's not just summing the real component of the basis.
For a complex input, we're doing all the same summation, and can visualize it the same, except $(AD + BC)$ no longer ends up being zero.


float image to integer type image I have a multispectral satellite image and I would like to "transform" each of the image bands into one with integer values. For example to stretch? the image between 0 and 255. Not sure how to proceed with this. I also thought I could make a histogram and replace each pixel value with the bin edge in which it falls into. I'm not sure. I'm currently working in python. Any thoughts on this?

You should check the Envi documentation (a standard tool for working with multi/hyper spectral images) on the stretches it uses to fit multi/hyperspectral data into 8 bits for display.

Here is the page on some of the stretches used, which apply an affine function to the pixel values and clip the values to the lowest and highest displayable value. In particular, look at the optimized linear stretch, which is the default one Envi uses to display things.

I think your description make it harder to understand the question.
It would be easier if you show the range of float image you have.
Is it between 0 and 1?

if it is, so you want to cast it into 0-255 integer. right?
just multiply your image by 255 and round it to integer.

Note that, this casting will make your image worse (errors from rounding)