OFDM RF signal waveform From the Wikipedia article:
"An inverse FFT is computed on each set of symbols, giving a set of complex time-domain samples. These samples are then quadrature-mixed to passband in the standard way. The real and imaginary components are first converted to the analogue domain using digital-to-analogue converters (DACs); the analogue signals are then used to modulate cosine and sine waves at the carrier frequency,  f{c}, respectively. These signals are then summed to give the transmission signal, s(t)."

So, information is added to a pure sine (or cosine) wave. This causes the wave to have not just one frequency, but a range of frequencies within a bandwidth that can be seen in the frequency domain of the signal.
 Having found out about this only yesterday, there are a number of questions.
 I've added an image of a regular time domain OFDM signal found somewhere on the internet.

Does the modulated signal really pass through the very narrow (ideally, 20 MHz in regular LTE) bandpass filter at the receiver? At a center frequency of 2 GHz that would mean to bandpass only what is within 1.99 GHz and 2.01 GHz, which sounds impressive to me.
Wouldn't even the slightest distortion from the emitter move the frequency of the signal out of the passband?
At 2 GHz operating frequency (analog time domain signal, at the LNA), there are 2*10^9 cycles in a second. If  the number of carriers is 2048, and a symbol is a string of 2048 cycles or slightly more including the cyclic prefix, doesn't it mean that the theoretical maximum symbol rate is (2*10^-9)/(2048+prefix length)? This would mean that the signal is transmitted continuously.
This is not true, however, since the symbol rate is a function of the sampling time, Ts, which is 32.552 ns. So, I don't really understand what FFT samples has to do with the 2 GHz analog signal that is present at the receiving LNA. My reasoning is that I could digitize the signal with any number of samples, get a plot for each 2048-cycle sequence then apply FFT.

If the graph below is the LNA signal of 2 GHz LTE, what would the period of a cycle be in nanoseconds? There are about 20 cycles from what I understand. How many cycles of this type would a symbol take?

Later edit:
Regarding the third question, I've found an answer.
Digital processing of the signal is done in the baseband, which is independent of carrier frequency. It needs to be converted into the passband. These links have more information:
IFFT and OFDM upconversion

baseband and passband modulation

Later edit 2: it seems that there are many ways to convert the baseband OFDM into the passband. Constant Envelope OFDM is one interesting way to do it, the baseband OFDM signal is phase modulating the carrier. This is in contrast to regular OFDM where the baseband is somesort of amplitude modulating the carrier.
More info at "Constant Envelope OFDM Phase Modulation", for example.

(1) Yes, if the technology permits it for the power handling requirements of the design, the strategy is to filter as much as possible immediately at the front-end of the receiver prior to the LNA. The primary reason is jamming from stronger signals. Depending on receiver topology used we would be concerned with 2nd and 3rd order intermodulation products and receiver phase noise that can be induced on out of band signals if significantly stronger all creating noise in band that would limit our receiver sensitivity.

FBAR technology is an excellent choice for front-end filters combining low loss and exceptional rejection. I have used FBAR filters in the design of femtocell cellular base-stations where the constraint was to receive the entire band (in this case it was a duplexer so would be used to receive 1850-1910 MHz and transmit 1930-1990 MHz --- notice the 20 MHz spacing between these two bands that this technology was able to provide exceptional isolation up at 2 GHz). We use the low loss band-selecting filter due to linearity/intermodulation concerns with the LNA, however with FBAR one must also be very careful with inter-modulation effects especially when using these filters in multi-carrier designs. Here is a good online article with more details about FBAR and other filtering considerations.

http://www.linleygroup.com/mpr/article.php?id=11440

(2) The distortion you are referring to is concern with frequency offset. Frequency offset will always exist since the typical transmitter/receiver set do not use a shared clock. The two dominant sources of frequency offset are the offset of the clocks themselves and additional offsets under dynamic conditions due to Doppler.  The relevant standards will specify the maximum allowable clock offset for the design, but < 5 KHz is a reasonable maximum possible frequency offset at 2 GHz for lower cost technologies (this is driven by the ppm specification for the master oscillator used).  If I calculated correctly the Doppler at this frequency when the transmitter and receiver are moving apart from each other at 60 mph would only be 179 Hz. The combination of both of these numbers will be well within the passband of whatever front-end filter is used, as no such filter would have a brick-wall response.

(3) I think you are confusing the carrier (2 GHz) with the bandwidth at that carrier (20 MHz). The actual carrier frequency is irrelevant to how much data we can transmit (without consideration to path loss effects on signal power for the SNR component of capacity).  The theoretical maximum symbol rate is given by the Shannon-Hartley equation for maximum channel capacity which accounts for both SNR and BW:

$$C = B \log_2(1+ SNR)$$

Where:

C is the maximum capacity in B/sec
B is the bandwidth of the channel is Hz
SNR is the signal to noise power ratio (not in dB but ratio of power)

Note MIMO technologies exceed this by utilizing multipath to send multiple "channels" simultaneously.

For more info on Shannon-Hartley equation see:

https://www.gaussianwaves.com/2008/04/channel-capacity/


Confusion regarding a library of Digital Signal Controller I'll be using a NXP DSC with 56800Ex core and my aim is to implement a Elliptic IIR Filter in it.

The User Reference Manual of the 'General Digital Filters Library' (which is available here) has functions to implement IIR Filters.

But on page 3-19 shown above, it says

  The filter coefficients are calculated using Butterworth approximation.

Does this mean that we can only use Butterworth Filter? Also, as per my knowledge, we have to give the coefficients and therefore other filters can also be implemented.
So what exactly does the function do?
Can I implement a Elliptic filter using this function or not? If not, then what should I do use Elliptic filter?

The function you are refering to only implements an order-1 filter. There is another function for order-2 filter.

If you want to implement an order-3 filter, you will have to split it in an order-1 filter and an order-2 filter. There are many ways to do it, but cascading the filters, as we would do in the analog world, is a good solution, albeit not optimal IIRC.

As to your original question, since the function refers to an order-1 IIR filter, an order-1 Elliptic or Butterworth, will have the same coefficients if you design them with the same cut-off frequency.

However, starting with order 2 or more, the coefficients will not be the same.


Question on 2-D Convolution This was a question on one of my past year paper questions.

  Question:  Prove

  $$f(m,n)*g(m-r,n-s) = f(m-r,n-s) *g(m,n)$$

  where $f(m,n)$ and $g(m,n)$ are 2-D discrete functions, $r$ and $s$
  are integers and $*$ is the 2-D convolution operator.

Not sure how to start on this and need some guidance on it.

In addition to the above comment:

Show that the 2D Z-transforms of both sides are equal.

Step-by-step: If $F(z_1, z_2)$ is the 2D Z-transform of $f(m,n)$, then what is the 2D Z-transform of

$f(m-r, n)$?
$f(m-r, n-s)$?
$f(m, n) * g(m, n)$?
$f(m-r, n-s) * g(m, n)$?

(Yes, I realize that the 2D Z-transform may not be taught before this exam was presented, but it is a valid solution, nonetheless.)


Hilbert transform of a bandpass signal I have a bandpass signal x. I am taking hilbert transform of it. I want to reconstruct amplitude and phase of the signal.Therefore, i am using Hilbert transform.

hil = hilbert(x);

f = ?

amp = abs(hil);

figure, plot(f,abs(hil);

phase = angle(hil);

figure, plot(f,angle(hil)

Here how to define the frequency axis, so i can get the amplitude graph 's length equal to my bandwidth. Moreover, is this the correct way for amplitude and phase reconstruction? My bandwidth is here 40 KHz. Thank you everyone in advance

You need to use a bank of bandpass filters before applying the Hilbert transform to get per frequency band amplitude and phase change over time.


Does FFT magnitude of recorded sound depends on the bit depth used? Based on some reading, here's my understanding of how sound recording works:

Microphone outputs a voltage based on sound pressure. The microphone voltage is amplified and the amplified voltage goes into an ADC. The amplifier gain is such that the maximum microphone output is matched to the maximum input voltage of the ADC. The ADC converts the input voltage into a number, and this number's magnitude can depend on the number of bits. So let's say the max input voltage is 5V, this would be 255 for 8 bit, or 65535 for 16 bit.

When I record an audio stream using something like PyAudio on Python, the audio is composed of list of numbers that the ADC outputs, and that can vary based on if I used 8 bit or 16 bit encoding. Therefore, the FFT magnitude will vary too.

Please correct me if I'm wrong.

With very few exceptions any audio processing these days is done in floating point, i.e. the data that comes through the ADC is immediately converted to floating point numbers and normalized so that the max voltage corresponds to a value of 1.0. This eliminates any dependency of the scaling to the number of bits.

Fixed point processing is quite difficult and typically requires specialized data types and intensive scaling and noise management.

The output of an 8-bit ADC could be interpreted as the numbers $-128,-127, ...+127$ but it's almost always better to think of it as $-1, -127/128, ..  -1/128 , 0 ,+1/128 ... +127/128$


How is it determined that an adaptive filter has converged? Is it when the cost function such as the error signal or Mean Square Error (MSE) signal is minimised?

The adaptive filter is converged when the error is what they call Wide Sense Stationary, meaning the mean and variance of the error are unchanging over long time intervals.


FFT shows systematic deviations from analytical result I am calculating the numerical Fourier transform of an exponential decay exp(-|t|) and compare it to the analytically calculated result, a Lorentzian. I find that the numerically calculated spectrum contains systematically larger amplitudes than the analytical one and that this deviation increases with frequency. I'm using the python library numpy.fft.
My time array contains [-100., -99.99, ..., 0.00, ..., 100.] and the signal is correspondingly
[3.72007598e-44, ..., 1.0, ..., 3.72007598e-44]. I've plotted this "original" exponential together with the back Fourier transforms below.

I then calculate the FFT (and adjust the phase to get a purely real result). When I compare the numerical result to the expected Lorentzian
$2 \cdot \frac{1}{1^2 + (2 \pi f)^2}$, I get very good agreement at low frequencies that however gets worse towards high frequencies. The numerical result (absolute value = real part) is systematically larger than the analytical one (labelled "calc" in this plot).

This deviation is of course also visible when plotting the ratio and difference of these two curves ("calc" refers to the analytic result). The ratio of the quantities seems to be independent on the step width and the maximum of the time array.

The inverse FFT of the numerical result (shown in the first plot, labelled "ift") seems to agree even a bit better with the original function than the FFT of the analytically calculated Lorentzian (labelled "icalc"). This depends however on the choice of the time array; it can also happen that the inverse FT of the Lorentzian is slightly longer close to the exponential function.
My objective is to evaluate experimental data in the frequency domain by models in the time domain. I would like to replace the usually used analytical FT (which restricts us to very idealized models) by a numerical FT. The first step is of course to show that the different FTs yield the same result... Of course, the experiment will not measure the spectrum with infinite precision, but I would very much appreciate any help to get at least 1-2 orders of magnitude better in the agreement numerical / analytical.

You have aliasing
Your continuous spectrum is not band-limited, so you are violating the sampling theorem. Sampling in time results in periodic repetition in the frequency domain so your first mirror spectra are leaking into your base band spectra.


Will the capacity of a channel becomes unbounded if i increase its signal-to-noise ratio $S/N$ without limit? According to the Shannon-Hartley theorem the capacity $C$ of a channel which has a signal-to-noise ratio of $S/N$ and a bandwidth $B$ can be calculated to be $C = B \log_2 \left( 1 + \frac{S}{N}\right)$.
Heren if $B→\infty$, the capacity does not become infinite since,
with an increase in bandwidth,the noise power also increases. If the noise
power spectral density is $N_0/2$,then the total noise  powr is $N=N_0B$,so the Shannon-Hartley law becomes
\begin{align}
C &= B \; \log_2 \left( 1 + \frac{S}{N_0B}\right) \\
&=\frac{S}{N_0}\left(\frac{N_oB}{S}\right)\log_2\left( 1 + \frac{S}{N_0B}\right) \\
&=\frac{S}{N_0}\log_2\left( 1 + \frac{S}{N_0B}\right)^{(\frac{N_0B}{S})}.
\end{align}
Now
$$\lim_{x\to0}(1+x)^{1/x}=e$$
So now it becomes
$$C_∞=\lim_{B\to\infty} C=\frac{S}{N_0}\log_2e=1.44\frac{S}{N_0}$$
so here channel capacity does not become infinite; that means, it's bounded.

But what happens if we increase the signal to noise ratio without bound? Will that give unbounded capacity? Is it possible?

After searching I got this but don't know how is this is possible; can some one justify this?

By the sampling theorem, uniform samples taken at more than twice the bandwidth of a band-limited signal can be used to reconstruct the signal perfectly by sinc interpolation. For each different set of sample values there is a different band-limited function. In absence of noise, by increasing the bit depth of the samples one can store an unlimited amount of information per sample into the band-limited signal, and get it back simply by sampling and digitizing it again.


Meaning of tight frame in the context of texture synthesis I am currently reading "A parametric texture model based on joint statistics of complex wavelet coefficients", by Portilla and Simoncelli. In this article, they use the "steerable pyramid" to represent the image. This uses a set of filters to decompose the image in a set of something around 16 subbands in the frequency domain, varying in spatial frequency and orientation.

Problem: They mention (page 55) that their decomposition forms a "tight frame". I know the definition of a tight frame in general linear algebra setting, or say orthogonal wavelet decomposition. However, in those cases the set of vectors that form the frame has cardinality of at least the dimension of the space considered. It seems highly unlikely that the space of images they want to consider is at most 16, and thus highly unlikely that this is what they mean. Can anybody tell me what is meant by the notion of a "tight frame" in this case?

Background: I am a mathematician that recently started working in neuroscience and now reading the mentioned article, published in the international journal of computer vision. Therefore the problem is probably just a matter of different languages used in different fields.

Thanks a lot!

Since visoft requested i post the solution, here goes. It has been a while so not all the details are still fresh so feel free to ask anything if things are unclear!

Let's say we have an image $I$ and apply a filter, which is given by convolution by some $F$. The pixel values of $I$ convolved with $F$ are all the output of the inner product of $I$ with a translated (and less importantly: mirrored) version of $F$.

What i missed at the time was that each pixel should be seen as a coefficient of the new representation of the image. So instead of the 16 subbands i mentioned, the number of coefficients of our new representation is 16 times the number of pixels in the image.

Some notation: for pixel number $(a,b)$, let us write $F_{(a,b)}$ for $F$ translated $(a,b)$ steps (and mirrored) and $(I * F)(a,b)$ for the $(a,b)$-th pixel of $I$ convolved with $F$.

Also, let us call all the translated versions of $F$ our basisfunctions. One can check from the definition of convolution that if we convolve our filter output again with $F$, this results in multiplying every basisfunction $F_{(a,b)}$ by the coefficient $(I * F)(a,b)$ and adding them up. This results in
\begin{align*}
I * F * F &= \sum_{(a,b)} (I * F)(a,b) \cdot F_{(a,b)}\\
&= \sum_{(a,b)} \langle I, F_{(a,b)}\rangle \cdot F_{(a,b)}
\end{align*}
where the last equality follows by definition of convolution.

From this, one easily recognizes something that comes up using an orthogonal basis. Now the collection of $F_{(a,b)}$ for all $(a,b)$ is called a tight frame if the above mathematical expression is equal to $I$ itself! It is thus very similar to an orthogonal basis, except that one does not need the set $\{F_{(a,b)}\}_{(a,b)}$ to be linearly independent (if you are very sharp you will note that in this particular case it most likely will be linearly independent, but see the following).

In the paper, actually about 16 of these $F's$ are used. Let's number them $F_1$ up to $F_n$. This is called a tight frame if the following holds:
$$
I = \sum_n \sum_{(a,b)} \langle I, F_{n,(a,b)} \rangle \cdot F_{n,(a,b)}
$$

This expression is equal to (check this):
$$
I = \sum_n (I * F_n * F_n).
$$

Since convolution is equivalent to multiplication in the fourier domain, one can check this by going to the fourier domain, squaring all the $F_n$ pointwise and adding them up. That's simple to do. If that adds to an image that is equal to 1 in every pixel, we are dealing with a tight frame.

So concluding, i missed two things: convolving with one $F$ gives as much coefficients as there are pixels in the image, and convolving twice actually multiplies the basisfunctions with the coefficients and adds them up.

As I said, I wrote this quickly and it has been a while, so if anything is unclear feel free to ask again!


Derivation of Nyquist Frequency and Sampling Theorem I have been looking through different sites and questions over the internet about Sampling theory, but couldn’t find the clear definition of how nyquist frequency condition is derived? It would be great if someone could direct me to the derivation of the condition that sampling frequency should >= to maximum frequency in the signal

It really boils down to aliasing.  In continuous-time, if you have any two signals $x_1(t) = \sin(2 \pi F_1 t)$ and $x_2(t) = \sin(2 \pi F_2 t)$, then as long as $F_1$ and $F_2$ are distinct, the signals are, too.

But consider sampling at some time interval $T_s$, so that the sampled signals are $x_1(k) = \sin(2 \pi F_1 T_s k)$ and $x_2(k) = \sin(2 \pi F_2 T_s k)$.  If you have a condition where $2 \pi F_1 = 2 \pi F_2 \pm 2\pi n$ for any integer $n$, then the signals are not clearly not distinguishable, by the simple trigonometric identity $\sin \theta = \sin \theta + 2\pi$.  Do the math, and you find that the frequencies cannot be separated by more than the sampling rate to be distinguishable.

It's worse, though, because you also have the trigonometric identity $\sin \theta = -\sin \theta + \pi$.  When you consider that the only thing that allows you to distinguish one signal from another is the frequency, that means that if $F_2$ is a mere $\frac{1}{2 T_s}$ away from $F_1$, then $x_1(t)$ cannot be distinguished from $x_2(t)$.

This phenomenon of two distinct signals appearing as one after sampling is called aliasing.  It's pretty much the basis of the Shannon-Nyquist sampling theorem, and the Nyquist rate.  Because if you're sampling at $\frac{1}{T_s}$ and trying to distinguish signals further apart than $\frac{1}{2T_s}$, you just can't.

Hence, the Nyquist rate.


narrow band frequency response creation and general FR question i need some guidance with understanding frequency response characteristics.

i have the a signal A with is played on speakers and recorded via microphone which is signal B.
then i apply a window function, shift both signals to frequency domain and analyze the magnitudes,
but the peak for the most dominant frequency is not at nearly the same index. it jumps.

to do further analyzation of the signal my plan was to create the frequency response of my audio system, multiply my frequency domain signal A with the FR and the peaks should align.
am i roughly right with this?

if yes, additional question:
as my processing power is limited i want to analyze only a certain frequency band.
what would be the best approach to generate an frequency response for a certain frequency band? and if i consider only a few frequencies, is it possible, that other frequencies outside my "scope" affect my observed ones and make my results unpredictable again?

thanks for any help in advance!
cheers

Let's assume that speaker, free space and microphone can be modeled as linear systems and that the overall system transfer function is given by $H(\omega)$. Then the received signal is given by
$$
B(\omega) = A(\omega)H(\omega),
$$
where $A(\omega)$ is the signal sent to the speaker. So, yes, if you know the transfer function $H(\omega)$ you can calculate $B(\omega)$ and it should align with what you record with the microphone.

To determine $H(\omega)$ in practice you will have use discrete-time signals and play a unit impulse $a(n) = \delta(n)$, record the received signal $b(n)$ and calculate the transfer function with help of the discrete Fourier transform (DFT):
$$
B(k) = DFT[b(n)]
$$
Now the maximum frequency $\omega_m$ you are able to analyze depends on the sampling rate $\omega_s$ with which the samples of $b(t)$ are acquired after the microphone to obtain $b(n)$: $\omega_m = \omega_s/2$. Frequencies greater than $\omega_m$ can neither be measured nor calculated and will influence your result because of the Aliasing effect unless they're filtered with a low pass filter before sampling. An alternative approach is to transmit a signal $A(k)$ that only contains frequencies in the range you're interested. The transfer function for this limited frequency range can then be calculated by
$$
H(k) = \frac{B(k)}{A(k)}
$$
In this case frequencies lying outside this range won't affect your result as long as it's a linear system.

A comment on processing power: the sampling rate influences the clock rate in a real-time system and doesn't affect complexity otherwise. The cost of the DFT operation (implemented by the FFT algorithm) increases with the length $N$ of the recorded signal $b(n)$. The greater $N$ the greater the frequency resolution of $B(k)$.


Designing artificial signals I'm trying to design a sample signal that has two "components". I want the signal to look like the upper left figure, with the two components shown on the right. The time frequency spectrum is shown in the bottom left.

I'm using MATLAB to generate this signal, as I am running MATLAB code that uses these signals. I have been trying my heart out with a combination of chirps, sincs and cos functions, but I just can't for the life of me reproduce this sort of a spectrogram. What I have produced looks like this:

It's quite close I feel, but there is something missing. Can someone point me to a way to approach this?

EDIT:

This is the code I use:

fs = 5000;
x0 = 1;
t = -x0:1/fs:x0;
yi = cos(t).*sinc(t).*chirp(t,150,0.2,250,'quadratic');
yj = cos(t).*sinc(t).*chirp(t,50,0.4,150,'quadratic');
subplot(2,2,1)
plot(t,yi+yj)
subplot(2,2,2)
plot(t,yj)
subplot(2,2,3)
spectrogram(yi+yj,512,511,256,fs,'yaxis')
subplot(2,2,4)
plot(t,yi)

In the code below I have manually created chirps where the frequency does not increase linearly, but behaves in a cosine-shape, as suggested by Atul Ingle . This waveform resembles the shape of the lines in the spectrogram nicely.

fs = 5000;
x0 = 1;
t = -x0:1/fs:x0;
yi = cos(t).*sinc(t/1.5).*cos(2*pi*500*(1.8-cos(2.2*t)).*t);
yj = cos(t).*sinc(t/1.5).*cos(2*pi*500*(0.85-0.5*cos(2.3*t)).*t);
subplot(2,2,1)
plot(t,yi+yj)
subplot(2,2,2)
plot(t,yj)
subplot(2,2,3)
spectrogram(yi+yj,512,511,256,fs,'yaxis')
subplot(2,2,4)
plot(t,yi)


given a mixture of two signals, how to separate it to its base components? I have a signal that contains two components. How can I separate it two its two components in matlab?
Please note that the location and shape of components are different in each experiment and we do a lot of experiments. It should automatically detect the two components based on the frequency contents. I have added this comment because some people suggested
c1 = x(1:80);c2 = x(81:end);

I do not mean such separation.

I would suggest you take a look at 'independent component analysis', and 'blind source separation' tasks. In general, given a signal of mixtures of components to be separated to its contents well addressed by optimization in complicated cases. Here you show one example, but probably in each experiment you have something very different. I will also take a look at my files to see if still I have relevant codes to solve it or not.


Proof of a CTFT property How do I prove that the absolute value of CTFT of a positive continuous time signal is less than or equal to its value at central frequency(0)?

If $f(t)$ is a real non-negative function with Fourier transform

$$F(\omega)=\int_{-\infty}^{\infty}f(t)e^{-j\omega t}dt\tag{1}$$

the following must hold:

$$|F(\omega)|=\left|\int_{-\infty}^{\infty}f(t)e^{-j\omega t}dt\right|\le\int_{-\infty}^{\infty}|f(t)|dt=\int_{-\infty}^{\infty}f(t)dt=F(0)$$

since $f(t)=|f(t)|$.


Approximation and IIR Filters Consider the discrete-time time-invariant system with input $x[n]$ and
output $y[n]$ satisfying
$$y[n] = \sum_{k=1}^5{x[n-k]}$$
Consider approximating the desired system with a second-order IIR
system with system function
$$H^{'}(z)= \frac{1}{1+a_1z^{-1}+a_2z^{-2}}$$
Use the following error criterion:
where $h_d$ is the desired impulse response.
$$E = \sum_{n=-\infty}^{\infty}\left\lvert h_d[n]+a_1h_d[n-1]+a_2h_d[n-2]\right\rvert^2$$

How can the particular error function be useful in solving the system of equations or the desired impulse response?
Do I need to take derivative w.r.t both $a_1$ and $a_2$ and make it 0 to get the system of equations?
Now this is what I am getting:
$$0 = \sum_{n=-\infty}^{\infty}( h_d[n]h_d[n-1]+a_1h_d[n-1]h_d[n-1]+a_2h_d[n-2]h_d[n-1])$$
and
$$0 = \sum_{n=-\infty}^{\infty}( h_d[n]h_d[n-2]+a_1h_d[n-1]h_d[n-2]+a_2h_d[n-2]h_d[n-2])$$
Now how can I solve these system of two equations and how can I approximate the desired impulse response from that?

I think it's instructive to try to understand how this error function is derived and why it makes sense. First, the desired impulse response $h_d[n]$ is given implicitly by the input-output relation of the desired system:
$$y[n]=\sum_{k=1}^{5}x[n-k]=\sum_{k=-\infty}^{\infty}h_d[k]x[n-k]\tag{1}$$
The given error function is the error function minimized by Prony's method for the design of IIR filters. We try to approximate a given transfer function $H_d(z)$ by an IIR filter $H(z)=B(z)/A(z)$:
$$B(z)\stackrel{!}{=}H_d(z)A(z)\tag{2}$$
For the given example we get (with $B(z)=1$)
$$1\stackrel{!}{=}\left(h_d[0]+h_d[1]z^{-1}+\ldots+h_d[N-1]z^{-(N-1)}\right)\left(1+a_1z^{-1}+z_2z^{-2}\right)\tag{3}$$
For equation $(3)$ to be satisfied, we need all coefficients associated with negative powers of $z$ to vanish:
$$\begin{align}h_d[1]+h_d[0]a_1&\stackrel{!}{=}0\\
h_d[2]+h_d[1]a_1+h_d[0]a_2&\stackrel{!}{=}0\\
h_d[3]+h_d[2]a_1+h_d[1]a_2&\stackrel{!}{=}0\\\vdots\end{align}\tag{4}$$
In practice, we can solve $(4)$ in an approximate way by minimizing the sum of the squares of the left-hand side, leading to the given error function
$$E=\sum_{n=-\infty}^{\infty}\big(h_d[n]+h_d[n-1]a_1+h_d[n-2]a_2\big)^2\tag{5}$$
Finally, the optimal coefficients $a_1$ and $a_2$ are obtained by taking the derivative of $(5)$ w.r.t. $a_1$ and $a_2$ and equating them to zero. This results in two linear equations with two unknowns.


Pilot symbols in IEEE 802.11a In IEEE 802.11a standards, there are many fields. Of which is the Data field which has 52 available sub-carriers used for transmission of coded bits. If the subcarrier locations are numbered from 32,-31,...0,1,...31.

I have read in the book that

that 48 out of 52 subcarriers are used for transmission of data bits. While 4 pilots are used as pilots for phase and frequency tracking and training. The pilot subcarriers are located at -21,-7,7,21 with values 1 (at -21,-7,7) and -1 at 21.

I am just wondering, usually the preamble (fields before the data field) is used to estimate the channel by sending a pilot symbol (using ALL subcarriers). Then my questions are

1) what is the purpose of sending these 4 pilot tones?

2) I am guessing phase and frequency tracking is different from channel estimations?

3) My last question is do we send the pilot values on the predefined pilot subcarriers ALONG with the data bits? so basically is each OFDM symbol made up of pilot and data tones? This is different than the procedure in preamble where  ALL subcarriers are used to send one pilot symbol, right?

Thanks

A frequency offset $f_\mathrm{off}$ between the carrier frequency and the local oscillator (at the receiver) results in a linear phase, i.e. the received signal contains a factor $\mathrm{exp}(j2\pi f_\mathrm{off}/f_\mathrm{s}n)$, where $f_\mathrm{s}$ is the sampling frequency and $n$ is the discrete time. If the frequency offset is sufficiently low (e.g. because it has been estimated and compensated for in an earlier step) the linearly changing phase can be assumed constant during one OFDM symbol. This means that all samples of the received OFDM symbol have a constant phase offset that can be estimated using the pilot subcarriers. As the phase is changing slowly, this has to be done for every OFDM symbol individually. The pilot symbols can also be used to estimate phase noise.
Yes, see 1. Frequency and phase offsets are time varying effects while the channel is assumed to be constant for at least some OFDM symbols.
Yes, pilot tones are sent in every OFDM symbol, along with data subcarriers.


Feature detector vs dense detector What is the difference between a feature detector and a dense detector in computer vision? Are they the same or a subset of one another?
For the following algorithms (implemented in OpenCV), which of them are feature detectors and/or dense detectors?

SIFT
SURF
AGAST
KAZE
AKAZE
FAST
BRISK
ORB
GFTT
HarrisLaplace
StarDetector

For class matching, Hietanen in [29] compares several
binary descriptors and SIFT with different detectors, including
a dense grid. SIFT preformed better than other descriptors and
dense grid responds very well.
This result leads us to test dense detector for different
descriptors…

It seems that a dense detector is one that uses a dense grid on which to evaluate some feature metric.
This is opposed to "classical" feature detectors that search the image for points with features.
Frankly, spanning a grid over an image to reduce it to a lower resolution does sound like what I'd call a subsampling and image transform, but maybe I'm misunderstanding this!


How can I generate a sine and cosine function using an exponential function? I have been trying to generate sine and cosine sampled graphs by a given sampling frequency.

fc = 4092; % Carrier frequency Hz

fs = 16368; % Sampling freq Hz
ts = 1/fs; % Sampled time

nn = [0 :4091] % Array with numbers from 0 to fc

y1 = exp(i*2*pi*fc*ts*nn); % Exponential equation

cosine = real(y1);
sine = imag(y1);

subplot(2,2,1);
plot(cosine);
axis tight;
grid on;

subplot(2,2,2);
plot(sine);
axis tight;
grid on;

If you have defined a sampling frequency then you should use it in defining your time grid "nn", otherwise you are just incrementing in digital steps.

The grid "nn" is basically a digital grid rather than representing true sampling based on your defined sampling frequency. Also the range of your time grid should be in mutliples of the time period of sinusoid (it could be anything but if you want to define it in that terms then define it in terms of time period of sinusoid and not the frequencies itslef).I have modified your code below to view 5 cycles of the sinusoid. Hope it helps

fc = 4092; % carrier frequency Hz

fs = 16368; % sampling freq Hz

nn = [0 : (1/fs):5*(1/fc)] % array with numbers from 0 to fc

y1 = exp(i*2*pi*fc*nn); % exponential equation


Can FFT Substitute Correlation for This Case? I'm trying to determine frequency of a repeating pattern in a time series. The sampling rate 100Hz and the pattern repeats itself with periods between 0.5 second and 1.25 second. This corresponds to a frequency range between 0.8Hz and 2Hz. At each second, I want to find the period of the pattern with 10msec precision, i.e. I want to differentiate 0.7sec(1.4Hz) and 0.6sec(1.6Hz).

I can attain this goal using correlation: For possible windows lengths (periods) I calculate correlation between two successive windows with same length. The window length with highest correlation gives the period of pattern.

Can I obtain the same results by using FFT?

As far as I see if I use 128-point FFT with this 100Hz series, then it will give me 100/128=0.8Hz (1.25sec), 1.6 Hz (0.62sec), 2.4Hz (0.42 sec),... components which obviously won't meet my requirement of 10msec period precision.

I know I can increase precision by using 256-point FFT but this corresponds to 2.5-second window for 100Hz and I want to determine pattern frequency for each window of ~1second.

How did you plan to use the FFT in the first place?

If it is to naively look for peaks in the magnitude spectrum and deduce the frequency from it, this is a bad idea. As you have seen, this has a limited resolution in the lowest frequencies. Resolution limitations aside, looking for a peak in the FFT is a poor way of estimating periodicity since the true fundamental frequency might not be the one with the most salient amplitude.

However, you can use FFT as an intermediate step to efficiently compute the autocorrelation function (FFT with zero-padding of the input to avoid wrapping -> square of magnitude of coefficients -> IFFT). Looking for peaks in the autocorrelation function is a more reliable period estimator than FFT peaks.

Thus, given your requirements, you could use sliding windows of 256 samples, with a 61% overlap, compute the autocorrelation through the FFT, look for the autocorrelation peak in the [50, 125] samples range and use this as an estimate of the period. Because of the 61% overlap, this will give you one period value per second. Observe here that it is possible to use large analysis windows (256 samples, or more) and make them overlap so that we can still get a detection for every 100 sample. The downside is that if the period of the signal rapidly change, your detection will be "smudged".

I see a few other approaches for your problem that do not suffer from any resolution compromise on the temporal axis:

If the waveform of your input signal is known (square, sawtooth, sine...), you could try locking an oscillator of similar waveform to it through a PLL, and use the PLL frequency as your estimator.
If your signal is purely sinusoidal, you can fit a complex exponential to it (through eigenanalysis of the autocorrelation matrix).
An approach that would give a very high resolution (one decision per sample) would be to run your signal through a bank of 75 comb filters in parallel, with delays of 0.5, 0.51, 0.52 .. 1.25 seconds. Compute the energy of the output for each of these 75 signals ; the one with the highest energy gives you the period.

Note that your problem has very similar requirements to that of tempo detection (onset detection yields an envelope signal at 100 Hz or the likes, we aim to get a tempo map in the 60 BPM - 180 BPM, ie 0.33 - 1Hz range, with a temporal resolution of 1s). Historically, the comb filterbank approach was used ; while most recent solutions are based on autocorrelation or similar autocorrelation-based pitch estimators (yin) ; with analysis windows in the 5 - 10s range.


Method to generate binary sequences with desired cross-correlation and autocorrelation properties I am looking for a method to generate two binary sequences of length $N$ with the following properties:

Very good cross-correlation for all shifts.
Good enough autocorrelation for all shifts except zero delay (good enough means $\gg 1$ and $\ll N$).

I can even relax the constraint on the cross-correlation as follows:

Very good cross-correlation for $\pm M$ shifts around zero delay, and no matter for the other shifts.
Good enough autocorrelation for all shifts except zero delay.

Following publications of Dr Welch and Dr Sarwate, it is possible to obtain a set of $K$ sequences of length $N$ and norm $N$ with the following properties :
$$\max(θ_a,θ_c)≥N\sqrt{\frac{K}{NK-1}}  \approx \sqrt{N},$$
and
$$\theta_c^2 + \frac{N-1}{N(K-1)}\theta_a^2 \geq N,$$
where $\theta_a$ is the maximum autocorrelation and $\theta_c$ is the maximum crosscorrelation.
From that, I understand that it is possible to find two sequences whose cross-correlation is ideal ($\theta_c = 0$) and whose autocorrelation is good enough ($\theta_a \approx \sqrt{N}$).
However, I am not aware of a method or an algorithm to generate such sequences.
Could you provide some directions or point to a paper that would help?
Besides, by limiting the constraint on the cross-correlation to some shifts only, can we expect to reduce noticeably the maximum autocorrelation? Are there methods to generate such sequences?

Flattered though I am that a short note (D.V. Sarwate, "Bounds on crosscorrelation and autocorrelation of sequences", IEEE Trans. Information Theory, vol. IT-25, pp. 720-724, 1979) that I wrote over 40 years ago is still receiving some notice, the OP's question is way too broad for this forum; people (not me!) have written whole books on this subject.
Neither Lloyd Welch nor I ever claimed that there exist sequence sets that meet these bounds exactly or almost exactly. Rather, every set of $K$ sequences of period $N$ must be such that $(1)$ and $(2)$ in the OP's question are satisfied. In particular, the OP's belief

From this, I understand that it is possible to find two sequences whose cross-correlation is ideal ($\theta_c = 0$) and whose autocorrelation is good enough ($\theta_a \approx \sqrt{N}$).

is a false understanding. If two sequences are uncorrelated (that is, $\theta_c = 0$), then $\theta_a$ must be at least $\sqrt N$, but it is not the case that there must exist a pair of uncorrelated sequences with $\theta_a \approx \sqrt N$. My paper cited above gives examples of a set of $N$ sequences (the rows of a length-$N$ DFT matrix) for which $\theta_c = 0$ but $\theta_a = N$ which is much larger than the desired $\sqrt N$. Unfortunately, while smaller subsets do still satisfy $\theta_c = 0$, the maximum off-peak autocorrelation continues to cling to $\theta_a = N$. The paper also described the Frank-Zadoff-Chu sequences for which $\theta_c = \sqrt{N}$ and $\theta_a = 0$. Both sets are optimal with respect to these bounds.  The sequences in both sets are complex-valued, not binary.
For sets of $N$ binary sequences of period $N$, a bound due to Sidel'nikov tells us that $\theta_\max$ must be at least as large as $\approx \sqrt{2N-2}$. $\,\,$ Gold sequences of period $N = 2^{2n+1}-1$ essentially meet this bound ($\theta_\max $ equals $\sqrt{2(N+1)}+1)\approx \sqrt{2N-2}$, but for periods $N = 2^{2n}-1$, $\theta_\max$ is $\approx 2\sqrt{N}$ for Gold or Gold-like sequences. If $\theta_\max\approx \sqrt{N}$ is non-negotiable as a criterion for sequence design, then for $N=2^{2n}-1$, the small set of Kasami sequences has $\sqrt{N+1} = 2^n$ sequences with $\theta_\max= \sqrt{N+1}+1 = 2^n+1$. For more details, see D.V. Sarwate and M.B. Pursley, "Cross-correlation properties of pseudorandom and related sequences," Proc. IEEE, vol.68, pp.593-619, May 1980.
With regard to other queries in the OP's question, some work has been done towards the construction of of sequences whose autocorrelation values are small for small nonzero offsets (while larger bounds for large offsets are OK), but typically, small crosscorrelation values for small offsets while allowing larger values for larger offsets has not received much attention. The former problem is of interest because it makes fine-tuning of the synchronization process in spread-spectrum systems easier; the latter might be more of a niche problem for special systems with nearly synchronized times-of-arrival.


Discrete-time Fourier transform of $\frac{\cos(\frac{n\pi} 6)}{(n+3)\pi}$ Find the Discrete-time Fourier transform of $\frac{\cos(\frac{n\pi} 6)}{(n+3)\pi}$

I thought of making it to be a sinc, but at the bottom there is $n+3$ and if I replace $n+3$ then I don’t know how to get rid of $-\pi/2$ in the new $\cos((\pi/6)m-0.5\pi)$

Could someone through me a hint?


Bandwidth of time-varying harmonics of speech For speech signals, do the harmonics have a constant bandwidth above some threshold (arbitrarily, -6 dB  dropoffs)? Or does this bandwidth decrease with the upper harmonics? If this is not known, how can it be measured? (I could imagine one way would be find -6 dB dB points around peaks of a given frame of samples).

If the process that creates a spectrogram includes the use of a windowing function, then the artifacts of that windowing will have a nonzero bandwidth associated with that particular window function (including its length).

If a signal is modulated, then it will have a non-zero bandwidth.  Any tiny variations in the mouth cavity shape or in the amplitude of the glottal closures or in the time period between glottal closures of voiced speech could be a cause of modulation of the harmonics.  These factors seem unlikely to be a fixed constant for all humans at all times.  Which might lead one to doubt a hypothesis of fixed bandwidth.


Filtering method that minimizes ringing in time domain, implementable in Python I am acquiring data through an oscilloscope, which I am then filtering. I am doing all of this in Python.

I need a bandpass filter with stop bands from 0 to .4MHz, and 4MHz and up. My pass band needs to be in the remaining region and would ideally be as large as possible.

I checked out a text from the library to help design a filter. The first filter I built was a Butterworth filter. However, the time domain response was awful for any decent frequency range, as I need an order 5+ filter.

I am now using the Remez (in scipy.signal) algorithm to generate my Transfer function coefficients. It looks All right, but I'm wondering if there is a better filtering option for flat time domain a fast frequency response?

It would appear by your question that you are new to filter design, which makes it difficult to answer questions about your filtering needs, so you might want to try this link.  http://www.iowahills.com/

This site has free FIR and IIR filter design programs that make it point and click simple to design and compare filters. The programs make it easy to see the how the various filter parameters such as pole count, tap count, response type (Butterworth, Chebyshev, Adjustable Gauss, etc.), linear phase, minimum phase, etc., affect a filter's response in both the time and frequency domain.

These programs can't replace Python, (they only design filters and filter data on the program's test bench), but they will make it quite easy to choose a filter design to use within Python.


Error reading a .PCM file I want to convert a .wav file with a sampling frequency of 44100 Hz to a 16 depth .pcm.

I don't why I'm getting those peaks at the beginning of  .pcm plot (second figure below 2). If you could explain and tell me how I can correct it, I would very much appreciate it.

This is my attempt:

% determine the least common multiple (lcm) of fsin and fsout
fsin  = fs;
fsout = 22050;
m     = lcm(fsin, fsout);
% determine the up and down sampling rates
up   = m/fsin;
down = m/fsout;
% resample the input using the computed up/down rates
x_22 = resample(x, up, down);
audiowrite([a_filename(1:11),'_22050','.wav'], x_22, fsout);
precision = 'int16';
fidr      = fopen([a_filename(1:11), '_22050','.wav'], 'r');  % open .wav file to read
fidw      = fopen([a_filename(1:11), '_22050','.pcm'], 'wb'); % open .pcm file to write
w         = int16(fread(fidr, inf, precision));% read wav file
fwrite(fidw, w, precision);
fclose(fidr);
fclose(fidw);

fidr2     = fopen([a_filename(1:11),'_22050','.pcm'], 'r');
[data, ~] = fread(fidr2, 'short');
fclose(fidr2);

% plot
figure(1)
set(1, 'color', 'w')
subplot(311),plot(x)
grid on, box on, axis tight, title('.wav (f_s = 44100 Hz)')
subplot(312),plot(x_22)
grid on, box on, axis tight, title('.wav (f_s = 22050 Hz)')
subplot(313), plot(data)
grid on, box on,  axis tight, title('.pcm (16 bit depth)')

  These are the resuts:
  Here are all the steps of the conversion:

This is the zoom in on the thrid subplot to see the peaks:

That sure looks like your are including the .WAV file header as part of your data.

When you read your file in, and write it out, using the file operations, you are not accounting for the header.  The header can vary in size (it's actually three RIFF headers), but is usually 44 bytes long.

These are the headers in my own (C++) words:

struct FirstHeader
{
char     RiffID[4];
int      RiffLength;
char     WaveID[4];
char     FormatID[4];
int      FormatLength;
};

struct SecondHeader
{
short    Always0x01;
short    TrackCount;
int      SamplesPerSecond;
int      BytesPerSecond;
short    BytesPerSample;
short    BitsPerReading;
char     Filler[16];
};

struct ThirdHeader
{
char     DataID[4];
int      DataLength;
};

Do a search on "Wav file header" and you can get lots of details.  If you don't need any of the parameters, I suppose just looking for the "data" in the header and the next 4 bytes are the UINT32 value of the length.  So your data starts the byte after.

Otherwise, follow the RiffLength and FormatLength to get to the DataID field.

The file can be mono or stereo, 8 16 or 24 bit, etc.


time frequency localization using wavelet transform I am currently doing analysis on Photoplethysmograph (PPG) data and I want to know the frequency (heart rate) at every time point so a windowed FFT might not be the best option. I am looking at wavelet to generate frequency and time information. I have been working with matlab example code however I have trouble determining what is the best wavelet to use for this application. I do now know much about wavelet. Does the wavelet used being dictated by the shape of my signal. In my case PPG signal which looks like this:

Or some other consideration is needed when choosing wavelet? How critical is the choice of wavelet?

Thanks,
Kelvin

What information are you trying to extract from your signal?

  I want to know the frequency (heart rate) at every time point

If this is the information you want then any sort of frequency analysis is unlikely to be very useful. It will show you that you have a 1 or 2 hz periodic signal with a particular frequency profile, but you already know all of this information so that won't be especially illuminating.

You are probably better off with peak detection, or local minima/maxima detection to extrapolate beats.

If you can provide a sample of your data, and/or more details about the exact information you are trying to extract then it will help us to find a more precise solution to your problem.


Generating drum sample in programming language How would you generate the sample of a drum, say "tom". Using the programming language and it's functions like addition, multiplication, sine, sqrt etc. I tried to generate the sine wave and then sharply cut the volume and add some white noise. This trick gave me something somewhat drumish, but far from perfect.

There are different approaches to generate artificial (I guess that's what you're aiming at here) drum sounds. The approach that you describe is certainly a good way to start - i.e. look at an existing sample that you're trying to approximate and find the properties that make up it's main characteristics. For example

is it tonal (e.g. like a tom or bass drum) or noise-like (e.g. like a high-hat) or a mixture (e.g. like a snare)?
what's the frequency/spectrum of the tonal part
what's the bandwidth of the noise part
what's the attack / decay behaviour
etc.

Knowing these important aspects of many percussive sounds, you can try to modify - e.g. - a sine wave by changing the frequency and amplitude over time, adding (filtered) noise and so on. Some experimenting will yield sounds similar to that of the classic TR-808. This machine may also be a nice source of inspiration on how to approach drum modelling using basic synthesis elements.

Another - somewhat more involved - approach is called physical modelling. This method starts from trying to emulate and - as the name says - physically model the elements that are involved in actually generating a sound. In the case of drum synthesis this comprises, for example:

the drum head
the stick that is used to hit the head
the drum shell.

A thesis on exactly this topic can be found here.

Resources on how to actually program these models can for example be found in the communities of Pure Data, Csound or SuperCollider as these are programming environments that have been specifically designed for audio/music signal processing (for example, see here for a short tutorial on drum synthesis in Pure Data).

For some theoretical background, I have been recommended the musimathics books, but have not read them personally.


Advantages of Splines in Compressed sensing I am working in compressed sensing, basically in ECG reconstruction

Normally we were using biorthogonal wavelets for reconstruction.

Someone advised me to use exponential spline basis instead of biorthogonal wavelets and I tried it and the results were better compared to wavelets.

Even though I got better results I still don't know why it worked. I am going for presentation and I am expecting this question.

That's the old question:

You say it's better, and wonder why, and everyone can only ask you to mathematically define better.

As soon as you define a measure for processing quality, you will be able, typically with not too much linear algebra, derive based on the properties of your signal, your measure and your basis, why one is better than the other.


Calculating maximum signal frequency knowing its sampling rate? I was curious as to how I could calculate the maximum frequency of a signal knowing only its sampling rate (for example 125 nanoseconds). Does the sampling theorem have anything to do with it (sampling rate is at least two times the maximum frequency)? So is the max frequency half of the sampling rate?

The maximum possible frequency in the sampled signal is half the sampling rate, assuming that you are not bandpass sampling and/or there is an anti-aliasing filter. The question is whether the signal has content up to that rate or whether it stops at a lower frequency.

One way to determine this would be to take your sampled signal $x[n]$ and transform it to the frequency domain via FFT to get $X[k]$.

Then find the largest value of $k$ for which
$$
\left |X[k] \right| > t
$$
for $0 \le k \le \frac{K}{2}$ where $K$ is the length of $X$ and $t$ is some "threshold" below which you decide the signal level is effectively zero.


What is the relation between input and output PSDs given system transfer function $H(s)$ If I have the system transfer function $H(s)$ in the complex frequency domain, how would I relate the input/output power spectral densities?
I have come across the relation $P_{out}(f) = |H(f)|^2P_{in}(f)$ in the frequency domain, where $P_{out/in}(f)$ refer to the input and output PSDs. Would I be able to use this same relation in the complex frequency domain as $P_{out} = |H(i\omega)|^2P_{in}$? Although I suppose that would mean the PSD would be in complex frequency domain as well?
This is all very new to me so any clarification or resources that I could look at would be greatly appreciated.

If the system described by the transfer function $H(s)$ is stable, you can obtain its frequency response by substituting $s=j\omega$, and use the relation that you found:
$$S_Y(\omega)=S_X(\omega)\big|H(j\omega)\big|^2\tag{1}$$
where $S_X(\omega)$ and $S_Y(\omega)$ denote the power spectra of the system's input and its output, respectively.


How Does Mean Centering Affect the Result of Using SVD to Compress Images? I have been learning about using the Singular Value Decomposition to find low rank approximations to matrices.
I had an image which I converted to a matrix. I regarded each row of the matrix as a 'data point'.

I did two things:

Found the mean of the data points then subtracted the mean from each data point and calculated the SVD of the result
Directly calculated the SVD of the original data matrix

Then I calculated the low rank approximations of the results from 1. and 2. by replacing some of the singular values with zeros.
Finally I converted the results back to images.

Both images looked quite similar. They both looked like how SVD-compressed images should look (according to textbooks).

From a principal component analysis point of view, not mean-centering the data before finding the low-rank approximation shouldn't
work well. It is equivalent to calculating the principal components from the eigenvectors of $\frac{1}{n-1}X^TX$ (where $X$ is the data matrix) instead of the covariance matrix
of $X$. So why did it appear like mean-centering didn't matter? Or does it depend on the type of images?

I will try to illustrate why it is important to remove the Mean from data when doing the PCA (SVD is the tool to basically do the PCA approach to dimensionality reduction -> Compression).

If the data isn't centered you will have one (Most probably the most dominant) component (Direction) which will point to the center of the data.
Since all directions must be orthogonal, it means the next directions, which should be correlated to the data itself will have harsh constraint.

You want the directions to be chosen according to the data layout and not according to its position in space.


How to get the sum of unit steps and unit ramps from a discrete signal? I have the discrete-time signal depicted below.

This signal can be written as follows

$$y[n] = \{0, 0, 2, 4, 6, 6 ,6,0,0,\ldots\}\,,$$

or as the sum of shifted impulses as:

$$y[n] = 2 \delta [n-2] + 4 \delta [n-3] + 6 \delta [n-4] + 6 \delta [n-5] + 6 \delta [n-6]\,.$$

In this case, how can I obtain the sum of steps and ramps from this signal?

The answer is:

   $$y[n] = 2r[n-1] - 2r[n-4] - 6u[n–7]$$

To avoid mistakes, let us first divide the original sequences by $2$ (for the sake of linearity):

0 0 1 2 3 3 3 0 0 0 0$\ldots$

It can be interpreted as a piecewise linear sequence, so you can expect a solution with about  three degrees of freedom with ramps and steps. As these basic primitives are causal, you can use the method of deflation (or successive approximations). This is a quite common approach, used for instance in subspace decomposition, or in basis pursuit in signal processing. The idea is: find a first approximation that matches what you want approximately, subtract it from the data, and re-iterate the processing. The convergence is not always guaranteed, but it can give you insights or heuristics  for a deeper approach. Here, it is likely to work fine.

You can start from left to right, build the obvious, remove it from the data, and iterate the process.

The first 0 0 1 2 3 can be approximated with ramp $r(n-1)$. Subtract  $r(n-1)$ from the data, and get:

0 0 0 0 0 1 2 6 7 8 9$\ldots$

You can compensate the "novel" left-most positive  slope 0 1 2 with $-r(n-4)$. You now get:

0 0 0 0 0 0 0 3 3 3 3$\ldots$

and then the $-3u(n-7)$ finishes the job. Oh no, don't forget to multiply by two again.

Here, the exercise is relatively easy to manage. Such a morphological decomposition can be used on more complicated signals. If you want more details, you can check for instance Automatic decomposition of time series into step, ramp, and impulse primitives, Galati & Simaan, 2006:

  Time series data that can be modeled as linear combinations of
  weighted and shifted primitive functions such as ramps, steps and
  impulses are representative of many industrial, manufacturing, and
  business processes. Data of this type also are found in statistical
  process control, structural health monitoring, and other system
  diagnosis applications. Often, the existence of one or more of these
  primitive functions may be indicative of the occurrence of a specific
  process event, making their detection and interpretation of great
  interest. The human eye is an exceptional tool at this kind of pattern
  recognition. However, for processes that generate large amounts of
  data the human eye encounters difficulties related to speed and
  consistency necessitating an automated approach.


Moravec, Harris noisy window Harris and Stephens writes about the interest window of Moravec:
"The response is noisy because the window is binary and rectangular", and suggests applying a Gaussian window.
My Question: Why is the response not noisy after applying a Gaussian window? Is it because a Gaussian filter removes noise from a picture, or because a circular window is better for sampling somehow?
Thank you.

I've added a link to the paper you're quoting in the question.

The issue is that both the Moravec corner detector, and the change to it suggested first in that paper:

use the derivative of the image. ANY difference operation will tend to enhance higher frequency noise in the image.  The idea with the Gaussian smoother is to reduce this effect.  The two operations (Gaussian smoothing and difference operation) are often combined into a single operation: the difference of Gaussians.

As the paper suggests, the anisotropy of the Moravec detector is why the differential is used, rather than the up, down, and diagonal differences in the Moravec algorithm. The smoothing is introduced to counter the noise issue.


Convolution introduces NaNs I am performing convolution in MATLAB using the built in conv function. I am convolving an impulse response (length 500 samples)  with an 1800 sample long input signal. This input signal contains 450 NaNs at random indexes between sample 530 and 1285. The output signal of this convolution is length 2299, and now consists of 1255 NaNs.
I'm unsure why this happens, and how to avoid it?

I'm also not sure I've provided enough info, but if I can add anything else please let me know.

and how to avoid it?

Address the NaN in your input data. Don't "fix", "paper over it" or "replace by 0". Find the root cause for the NaN, understand what's happening and take meaningful corrective action.

NaN means your input data is bad. Doing anything with bad data is pointless since your output will be bad. Doing cosmetic adjustment just so the code doesn't blow up is equally pointless. You need to understand why the data is bad and address the root cause


Sine wave in frequency domain on signal analyzer The Fourier Transform of sine wave ($\sin(2 \pi A t)$) is given as :
$$\frac{1}{2i} [\delta(f-A)-\delta(f+A)]$$

This means that the Fourier Transform of the real function, $\sin(t)$ has an imaginary Fourier Transform (no real part).

How can we observe this Fourier Transform of sine wave on signal analyzer if it is imaginary with an "$i$" in it?

The thing is that I have a very high frequency signal (a few GHz) and so I am not able to see this signal clearly in time domain because the resolution of my signal analyzer is not that good. So I was trying to view it in frequency domain instead of time domain. Right now, in frequency domain I am able to view this signal as a single spike. So, should I conclude that the signal generated is a sine wave?

in frequency domain each spike/impulse represents to a sinusoid in time domain. so if in frequency domain the signal shows a single spike/impulse it means the input signal has a sinusoid in time domain.


Find the frequency response if i have the magnitude response? if i have the transfer function of magnitude response is there a method that i could calculate the frequency response?

For example the transfer function of the magnitude response is:

$ 3db \pm  3.5db  $ for $|ν|<0.1$

$ <-55db $        for  $|ν|<0.2$

The frequency response of a system can be represented in polar format, in which the magnitude and phase response are considered separately:

$$
H(\omega) = |H(\omega)| \angle H(\omega)
$$

With this representation, it should be clear that the magnitude response alone is not sufficient to characterize the full frequency response of a system; you have to know (or assume) what its phase response is also.


How can I remove random dots from an image? I recently received some scanned pictures which are in a pretty good state, with OK resolution, but they are randomly spattered with a bunch of black dots. These cover everything (clothes, faces, walls, windows, and so on) except for a few very bright areas. I also have no idea how they got there. Here is a sample:

From an information content perspective, I would expect that these dots could probably be removed as they likely have a quite distinct signature. (If they were evenly distributed, for example, I know that a spatial mask on their Fourier transform would do wonders.) Is there some algorithm to remove these artifacts?

OK, thanks to the fantastic Stack Exchange Ask a Question page, I have now a good idea of what they are and how to remove them. This is what I would describe as "pepper noise", as opposed to the traditional salt and peper noise. It can be removed by using a median filter, which essentially replaces each pixel by the median of the values of the pixels in some suitable neighbourhood around it. This is available directly in Photoshop, and as part of the Despeckle filter in the GIMP (see e.g. this resource). I hope this Q&A will be useful for other people with the same problem.


Beginner level: deconvolving Gaussian noise to find initial PDF noob question.

An initial process has an unknown pdf, P(X), which is then subjected to additive Gaussian noise.  Repeated sampling of the sum is performed so that

$y=x+noise$

is repeatedly measured and the pdf for y is determined, P(Y).  Is there anything that can be said analytically about P(X)?

I'd appreciate any simple references.

Thank you

You have a random variable $X$ whose pdf is $f_X(x)$. when another random variable $N\sim f_N(n)$ is added to $X$, then the pdf of the sum is just $f_Y=f_N*f_X$, where $*$ is convolution. This can be written in the Fourier domain as $F_Y=F_NF_X$. So if you know both $F_Y$ and $F_N$, and your $F_N$ is well-behaved, then $F_X=\frac{F_Y}{F_N}$.

After that, you can acquire $f_X$ by inverse Fourier transform.

Start by looking at the histogram of $y$, by which you can find $f_Y$. I think you can easily do the rest, given the above information.


Calculating transmitter power from spectral density I want to find out what is the real transmit power of my USRP 2922, unfortunately the only results I get when looking at the spectrum analyzer is the spectral power density.
These are the results of the my measurements (in LabView - transmitter parameters are: 1Mhz IQ Sampling rate, carrier: 2.4375GHz, 0 Tx gain, QPSK modulation; RX port on USRP receiver and TX port on USRP transmitter are connected with an attenuator):

Now, I'm wondering how to find the transmit power of the transmitter from this. The Gain control is for the RX Antenna gain, so it's not useful since it's not power of the USRP transmitter. Maybe I should integrate this value in the band? Alternatively for a sine of the signal, calculate the value from the amplitude? Any ideas?

You're right, to convert from a power density to a power, you integrate over frequency!
However:
Well, if neither of your transmitter nor receiver are calibrated, then you simply can't find the "real" transmitted or received power. The plots you have are "amplitude relative to maximum ADC input". You don't know what power "maximum ADC input" corresponds to without calibration.
End of story - there's not a single absolute power in your system that you can use as reference!
The solution here is to calibrate; i.e. send a signal, measure it with a calibrated spectrum analyzer, measure the same signal with your USRP, and then calculate the factor between these two powers. Due to linearity, that factor (==corrective additive term in dB) stays the same for all powers.


Phase Correlation - Poor Performance on Noisy/Blurred Images? I have successfully tested 1D phase correlation algorithm to determine vertical shift between two synthetic images.

When I moved to real images, however, it is not able to detect translation at all (the peak is located at 0, which is wrong result).

I have the following images:

And resulting phase correlation (Magnitude, Real, Imaginary):

The first scanline of the image is completely white, but the shift is obviously larger (20 pixels).

The expected result is white line on 20th row which happens only on synthetic images or light noise.

My algorithm is very simple - for each image column:

Compute 1D FT of source and target image columns (a=FT(A), b=FT(B))
Compute cross-power spectrum (cross_power = a *. conj(b) / |a *. conj(b)|) - *. denotes pointwise mutliplication, conj(x) denotes complex conjugate
Compute phase correlation (phase = IFT(cross_power))
Find maximum magnitude in every column of phase.
Find consensus peak location (e.g. median of detected peak locations)

Can you please advise me how to improve baseline phase correlation algorithm to deal with real-world (noisy) images?

Should I rather use NCC (Normalized Cross Correlation) instead of FFT-based phase correlation?

UPDATE

I was experimenting with zero padding to rule out errors introduced by circular shifting (only simple linear shifting of images is desirable) and tested this on original images from Wikipedia:

The single peak is clearly there as it should be:

However - if I perform slight smoothing (Gaussian blur) to reduce noise and actually improve the result, the phase correlation comes out totally mangled:

Here is the enhanced version - the original peak is weaker (why??) and there appeared new peaks around zero shifts (why??):

One-dimensional version

The one-dimensional version that you list won't work. When there is a large enough shift in images (more than one or two pixels in real-world images), there will be nothing relating the column pixels.

For an example of this, try:

I5 = rand(100,100)*255;
I6 = zeros(100,100);
I6(11:100,22:100) = I5(1:90,1:79);

So that we have I5:

and I6:

Then the one-dimensional phase correlation is just:

whereas the two-dimensional phase correlation is:

It's a bit hard to see, but there is a very high peak in the bottom-righthand corner of the image. No clear peak exists in the one dimensional version.

Why doesn't smoothing help? #1

What correlation is trying to do is to find "similar" variations in each image.  If the underlying signals are sufficiently random, then this will work well: the correlation of white noise with itself gives a really nice peak at the origin, and close-to-zero elsewhere.

Smoothing a "random" image with a Gaussian will have the effect of smoothing out the correlation your expecting --- spreading the energy in any peaks over a wider area.

Smoothing has the opposite effect of "pre-whitening" the image. Pre-whitening (as the name suggests) tries to make the image more like white noise --- which has the best form if we are doing correlation-based detection (in that the peak is well-localized).

What you are better off doing is using the matlab diff operation for a simplistic, but surprisingly effective, way of pre-whitening images.

See this example.

Why doesn't smoothing help? #2

Why does smoothing cause the extra peaks?

If you smooth each image with a kernel $k(x,y)$ then we get:
$$
h_a = g_a \star k\\
h_b = g_b \star k
$$
where $\star$ is convolution.

Now,
$$
H_a = K G_a\\
H_b = K G_b\\
R = \frac{H_a H_b^*}{|H_a H_b^*|} = \frac{|K|^2G_a G_b^*}{|K|^2|G_a G_b^*|}\\
=  \frac{G_a G_b^*}{|G_a G_b^*|}
$$

What I suspect is happening (though I am not sure) is that perhaps your kernel has near-zero values in the frequency domain, causing numerical problems?

If I apply a kernel:

K = one(5,5);

to my random images, then I get:

for the two-dimensional correlation, which makes the peak more spread out, but does not exhibit the problems you are seeing.


Transfer function block diagram Can any one help with the $y[n]$ and $x[n]$ relationship in this block diagram, I just keep have a $t[n]$ in my answer that I can't get rid off.

On my best try I got to $y[n] = 2t[n]-x[n-1]-y[n-2]+x[n]$

If you could also tell me how you got there it would be great.

Thanks anyone for help.

I prefer to solve such problems using the $\mathcal{Z}$-transform. First, write down the time domain equations:

$$\begin{align}\tag{1}
y[n]&=r[n]+x[n]\\
r[n]&=r[n-1]+t[n]+t[n-1]\\
t[n]&=t[n-1]+y[n-1]+y[n-2]
\end{align}$$

Taking the $\mathcal{Z}$-transform of these equations gives

$$\begin{align}\tag{2}
Y(z)&=R(z)+X(z)\\
R(z)&=z^{-1}R(z)+T(z)+z^{-1}T(z)\\
T(z)&=z^{-1}T(z)+z^{-1}Y(z)+z^{-2}Y(z)
\end{align}$$

After simplification you get

$$\begin{align}
R(z)&=T(z)\frac{1+z^{-1}}{1-z^{-1}}\\
T(z)&=Y(z)z^{-1}\frac{1+z^{-1}}{1-z^{-1}}
\end{align}$$

from which

$$R(z)=Y(z)z^{-1}\frac{(1+z^{-1})^2}{(1-z^{-1})^2}\tag{3}$$

follows.
Plugging (3) into the first equation of (2) gives

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

from which you finally obtain afters some algebra

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

In the time domain, Eq. (5) is equivalent to

$$y[n]-3y[n-1]-y[n-2]-y[n-3]=x[n]-2x[n-1]+x[n-2]\tag{6}$$


what is the difference between $X(j\omega)$ and $X(\omega)$ notation I am trying to understand Fourier Transform and Laplace Transform. What is the difference between  $X(j\omega)$ and $X(\omega)$ notation?

what is the meaning of $j\omega$ ? Is it represent frequency? If it is, what is the meaning of imaginary frequency?

Thanks in advance.

Both notations are common and correct. As pointed out by Yuri Nenakhov, the advantage of the argument $j\omega$ is that it coincides with the complex (Laplace transform) variable $s$ when its real-part is zero. Note that in the complex $s$-plane the frequency axis is the imaginary axis. So $j\omega$ has nothing to do with complex frequency (which makes no sense).

So if the Laplace transform $X(s)$ of a signal $x(t)$ exists, and if the imaginary axis is inside its region of convergence, then the Fourier transform is obtained by setting $s=j\omega$.

Note that this does not work in general! In general you can't get the Fourier transform by replacing $s$ with $j\omega$ and vice versa. Two conditions must be satisfied in order for this to lead to a correct result:

Both transforms must exist (in the sense that the corresponding signal $x(t)$ has a Laplace transform and a Fourier transform).
The imaginary axis $s=j\omega$ must be inside the region of convergence of the Laplace transform.

An example where replacing $s$ by $j\omega$ doesn't work, even though both transforms exist, is the step function:

$$\begin{align}&x(t)=u(t)\\\text{Laplace transform: }&X(s)=\frac{1}{s}\\
\text{Fourier transform: }&\hat{X}(j\omega)=\pi\delta(\omega)+\frac{1}{j\omega}\neq X(s)|_{s=j\omega}\end{align}$$


Detection of number of yarns in the fabric image I am a novice user in image processing and working on a project to determine number of yarns in a fabric image i.e. warps and wefts (Perpendicular and parallel yarns to x-axis).

I am trying to detect the yarns by taking the sum of columns and rows for warps and weft respectively but it doesn't seem working because the image is not very well focused and sharp. Also the method is much sensitive to the alignment of yarns in the fabric.

Please propose a better algorithm if possible for detection of yarns. Also if there is any possibility to make yarns straight for every captured image. I have tried IMROTATE but it creates the black regions in image which disturbs the sum of rows and algorithm gets failed.

Any help in this regard will be of much worth to me. Thanks.

A fourier transform will give you the period of any periodic features in the image - eg it will tell you there is structure event 3pixels horizontally and 5 pixels vertically, from this and the pixel scale you can calculate the number of yarns


Output of a linear time-invariant(LTI) system I am very confused about one of the questions I received during an exam. How do I solve this question?

This question could be rephrased as is $y(t) = T{x(t)} = x(t) ^ 2$ LTI system. This is a textbook case of a non-linear system. Here is a quick proof:
For a system to be linear the following must be true:
$T{a \cdot x(t)} = a \cdot T{x(t)}$ . i.e. if you input scaled signal into linear system you should get the same output as if you inputted that same signal but unscaled, and then scaled the resulting output. This is referred to as the scaling property.
We can use this test on the system in your exam question $y(t) = T{x(t)} = x(t) ^ 2$ :
$y(t) = T{a \cdot x(t)} = (a \cdot x(t)) ^ 2 = a^2 \cdot (x(t))^2$
$y(t) = a \cdot T{x(t)}= a \cdot (x_1(t)) ^ 2 = a \cdot (x_1(t))^2$
As we can see the results are not equal. Therefore the system is not linear, and so it is not LTI.
It appears that your teacher may have wanted to approach this problem in a different manner, so I apologize if this solution isn't helpful. However, this approach is completely valid and much quicker in my opinion. After getting more comfortable with it you can identify operators such as a square and immediately recognize that the system is not an LTI.


Nyquist Frequency on semi-unevenly sampled data I have a data set that has 'kind of' constant sampling rate - it switches between 1 min and 2 min. About 70% of the times, samples are taken every 1 minute, and about 30%, samples are taken every 2 minutes.

I tried using FFT, but it didn't work because samples aren't evenly spaced. Instead I applied Lomb-Scargle transformation, and obtained a good looking plot. Upper plot is a Lomb-Scargle periodogram, and lower plot is the original sample data with varying sampling rate (1 min or 2 min).

What should I use for Nyquist Frequency in this case?

the lomb scargle in matlab, (plomb) returns a frequency vector as the second output. i would be surprised if the routine you are using doesn’t do the same.  there is a common tendency for python signal processing libraries to be functionally equivalent to matlab conventions.

your plots show symmetry around a “center” frequency which is the same as what a FFT does for real inputs, except a FFT would not extend to start a third interval.  The high half frequencies(above $\approx$ .017 to .033)  are the negative frequencies in FFT speak. you can ignore them when you have a real signal.

There really isn’t a simple “Nyquist” frequency with nonuniform sampling.

The documentation of the library you are using should detail how frequency is treated. If it doesn’t, you can find a routine that does or do your own calibration by using test sines of known frequency.


Where to learn about "analog prototype filters"? Where to learn about "analog prototype filters"?

I've heard about them, but I'm unsure about what they really are and how they're constructed.

"Analog prototype" filters are well-known analog filters that have specific desirable properties. They include (but aren't limited to):

Butterworth filters: maximally-flat passband response
Chebyshev filters: sharper rolloff than Butterworth, equiripple in either the passband or stopband
Elliptic filters: even sharper rolloff, equiripple in passband and stopband
Bessel filters: maximally-flat group delay

These prototypes can be used to design digital filters that have approximately the same characteristics, for instance by using the bilinear transform.


Deconvolving the measurement system from the system response I'm preforming a system identification measurement, and would like to deconvolve the measurement system's response from the measured response.

The set up of the system is as follows; the signal is generated on a PC, the signal is then sent to the sound card, which in turn forwards the signal to the power amplifier, this drives a loudspeaker.  The loudspeaker propagates the sound energy which is captured by a microphone, and back into the sound card.  The recorded signal is then processed by the computer.

The signal from the sound card is also fed straight back to itself, which allows me to measure the response of the pc -> DAC -> ADC -> pc response.

I'd like to deconvolve this 'measurement chain' from the overall 'system chain'.  What methods could I use for this task?

Sorry for the long winded answer, but making a good acoustic measurement is hard. Here are some of the steps involved (and at some point, we'll actually get to the direct answer of your question).

First, you have to make sure that the system you are measuring is actual linear and time invariant (otherwise you can't deconvolve). This can be tricky on a PC if you use regular sound card drivers. These will route any signals through the Windows Kernel Mixer which frequently applies sample rate conversion and yanks buffers around in time in non- deterministic ways.I strongly recommend bypassing the Kernel mixer.

Second, you need to determine the length of the impulse of your system. A good estimate is the reverb time of the room that you are measuring. It's most convenient to choose a power of 2; for most rooms this will 16384 or 32768 at 44.1kHz or 48 kHz sample rate.

Third, create a periodic excitation of that length. This could either be a log sweep or (better) a pseudo random noise. The spectrum of the noise should be chosen so you get roughly constant signal to noise ratio in the frequency range of interest. That depends on the transfer function and on the background noise spectrum. If you don't know either of those yet, pink is a good start. Let's call one period of this signal x[n].

Fourth, wire up your system so that the left channel from your D/A goes into the loudspeaker AND into the left channel of your A/D. Wire the microphone into the right channel of the A/D.

Fifth, start the excitation (loop it or create a wave file with many repetitions of your noise signal). Carefully monitor all levels: make sure that the A/D sits at about 10 dB below clipping. Make sure the microphone pre-amp sits at about 10 dB below clipping. Make sure the power amp is not clipping and that the loudspeaker is not over driven.

Sixth, make sure that the room is as quiet as it can be. Close doors and windows. Shut off most things with a fan including any HVAC system. Evict all other people from the premises. If there are any ground loops, use isolation transformers and ground lifts as necessary. A good way of checking noise is to hook up a headphone amp to the mic output and listen to it through headsets. Any noise, hum, or other artifact that you can hear will also show up in the measurement.

Seventh, do the actual acquisition. Collect 12 periods with the excitation signal running. Visually inspect the results for something unusual (gaps, parts missing, drop outs, etc.)  Throw the first two periods away. Calculate the average over the other 10. Let's call the left channel y[n] (sound card) and the right channel m[n] (microphone).

Eighth, calculate the Fourier Transform of y[n]. This should be fairly flat with no zeros or areas of very low energy. This may not be the case as most sound cards have AC coupled inputs, i.e. there is some high pass filter and the value at DC might be very low.  Similarly there may also likely be an anti aliasing filter, so again you may have low energy or just noise at very high frequencies. If can you manually fix that (by adding a small amount of broadband energy), then do. If the spectrum is very flat in the frequency area of interest, you can simply replace it with a properly delayed unit impulse. If none of these will work, things get more complicated.

Ninth: deconvolute. Your transfer function can simply be calculated as $$H(\omega) = \frac{\Im \left \{ m(t) \right \}}{\Im \left \{ y(t) \right \}}$$
where $\Im \{ \}$ denotes the Fourier transform. This is direct spectral division without any zero padding or windowing. Since we are doing circular coherent acquisition we are also doing circular deconvolution.

The equation shows the problem discussed in step 8. We are dividing by a measured spectrum. Any zeros or areas of low energy or high noise in the divisor spectrum will result in garbage in the resulting transfer function. The results at the "good" frequencies will still be okay but you could not use it to, for example, calculate an impulse response. A simplified version would be to just divide by the excitation signal, i.e.
$$H(\omega) = \frac{\Im \left \{ m(t) \right \}}{\Im \left \{ x(t) \right \}}$$ and choose the excitation so that there are no zeros in the spectrum. However this will not allow you to deconvolve the D/A and A/D system response and it will also lose any absolute delay information. In many cases this is actually acceptable and many commercial measurement system do just that, so it depends on your application.

Tenth: verify that your measurement is good. Various tests should be performed:

Measure a few times and make sure that the result is the same.
Measure with the loudspeaker muted. This will give you a good estimate of the background noise spectrum. As a rule of thumb, you need at least 10 dB signal to noise ratio at all frequencies of interest.
Linearity test: measure with half the excitation gain and confirm that the resulting transfer function is the same.
A useful inline noise test is the following: Take your 10 periods and average 5 times over 2 periods each, then do a Fourier transform over 2*N samples. If your signal is noise free then all the odd bins should be zero. You can directly estimate the signal to noise ratio at any frequency as X{2*N+1}/X{2*N}


How do I optimize the window lengths in STFT? I have many EEG signals and I want to analyze them using linear methods such as STFT (Short Time Fourier Transform). In STFT , How can I optimize the analysis window length, to reflect the frequency spectrum of each analysis window in a proper way?

This is the classic "uncertainty principle" of the Fourier Transform. You can either have high resolution in time or high resolution in frequency but not both at the same time. The window lengths allow you to trade off between the two.

If you want to detect "events" in your EEG signal with a resolution of say 10ms, then this should be your window length. This will give you a frequency resolution of about 100 Hz.


Scaling property of Fourier Transform Problem 4.6(b) from Oppenheim, Wilsky & Nawab (2nd ed) reads:

Given that $x(t)$ has the Fourier transform $X(j\omega)$, express the Fourier transform of $x(3t - 6)$ in terms of $X(j\omega)$.

The way I proceeded,
$$ x(3t - 6) = x(3(t-2)) $$

$$ x(t) \longleftrightarrow X(j\omega)$$
$$ \implies x(t-2) \longleftrightarrow X(j\omega)e^{-2j\omega}$$
$$ \implies x(3(t-2))  \longleftrightarrow \frac{1}{3}X(\frac{j\omega}{3})e^{\frac{-2j\omega}{3}}$$
But the answer given at the back of the book is
$$  x(3(t-2))  \longleftrightarrow \frac{1}{3}X(\frac{j\omega}{3})e^{-2j\omega}$$

My question is, while applying the scaling property of Fourier Transform, why does scaling happen only to the $X(j\omega)$ term and not the exponential term, when even the latter forms part of the actual $X(j\omega)$ on RHS?

Ok you are mixing the order of the transforms. When you decompose the argument $3t-6$ into $3\cdot(t-2)$ , then you should first scale $x(t)$ by $3$ and then shift the scaled result by $2$, that's where you have chosen the wrong order. Let's show the correct order:

Let $x(t)$ be the original signal whose CTFT is $X(\omega)$, then define the following signals in two steps to get the final signal $z(t)= x(3t-6)$:

$$ y(t) = x(3t) \tag{1}$$ and
$$z(t) = y(t-2) \implies z(t) =  x(3(t-2)) = x(3t-6) \tag{2}$$

Notice how we replaced $y(t-2)$ by using its definition from the eq. (1).

Now applying CTFT to the above rows yields:

$$ Y(\omega) = \frac{1}{3} X(\omega/3) \tag{3} $$ and
$$ Z(\omega) = Y(\omega) e^{-j2\omega} \tag{4} $$

finally insert $Y(\omega)$ from eq. (3) into eq. (4) to get:

$$ \boxed{Z(\omega) = Y(\omega) e^{-j2\omega} = \frac{1}{3} X(\omega/3) e^{-j2\omega} }$$

By the way, you could have also used the alternate approach: first shift $x(t)$ by $6$: $$y(t) = x(t-6) \implies Y(\omega) = X(\omega)e^{-j6\omega}$$ and then scale $y(t)$ by $3$ to obtain $$z(t) = y(3t) = x(3t-6) \implies Z(\omega) = \frac{1}{3}Y(\omega/3) = \frac{1}{3}X(\omega/3)e^{-j2\omega}$$


How can convolution be a linear and invariant operation? I'm having a slight breakdown right now with a seemingly simple question. Say I have a system that convolves an input function with itself to produce an output function:
$g(x) = f(x) ∗ f(x)$
I've heard countless times that convolution is a linear operation and just assumed it was fact. However, I'm trying to understand how this is possible.
In my system above, it doesn't seem like this system could be linear since convolution contains multiplication. So multiplying an input signal with itself can't possibly be linear. However, it would appear to be shift invariant.
Am I correct in my logic? I seem to have a misunderstanding of convolution being a linear operation.

Your problem comes about because you're pondering the convolution of a signal by itself: $$g(x) = f(x) * f(x)$$
So, yes, if you had a block diagram where you collect f(x) from the æther and then somehow magically convolve it by itself (which can't be done physically because it would require time-reversing $f(x)$, which requires looking into the future), then yes, that would be nonlinear.
But that's not what convolution is about
Convolution is about taking an actual signal (probably what you mean by your $f(x)$), and running it through a system, who's behavior is defined by an impulse response.  An impulse response is, in a sense, a signal, but it's not an actual physical signal -- it's a full and complete description of how a linear, time-invariant system behaves in response to a signal.
Here's a fairly far-fetched analogy, that isn't coming out nearly as simple and clear as I'd like*: say you live on a farm, where goats weigh 1/3 as much as sheep, and you have a scale that reads out in goat-weights.  If you put three sheep on the scale, then the number you read is $3 * 3 = 9$.  What it looks like you have done is multiply three by three (a nonlinear operation if $3$ is a signal) and gotten 9.  What you have actually done is multiply a signal (three sheep-weights) by a fixed characteristic of a signal (one estimated-goat-weight / actual-goat-weight) and gotten an output.  Had you put four sheep on the scale, only one part of your calculation would change -- the "three goats per sheep" part would still be a three.
So in your original example, you're trying to convolve a signal "by itself".  That doesn't happen.  You may convolve a signal by an impulse response that happens to be equal to that signal -- and even I, in a classroom setting, may tell the class that I'm "convolving the signal by itself", only to have to make the above clarification.
You'll have occasion to find the autocorrelation of a signal, which looks a lot like convolution (it looks like $A_f(x) = f(x) * f(-x)$, in fact) and for which you sometimes trick a math package into doing because it doesn't have a built-in autocorrelation function.
You sometimes have occasion to find out what happens when you feed a filter with a signal that happens to match its impulse response -- but then you're feeding one distinct signal into a filter that's got eerily similar characteristics to the signal**.
You may even build a system that involves delay lines and multiplications and whatnot that "convolves a signal by itself" -- but then you're building a filter on the fly that gives itself an impulse response that matches a delayed version of a collected signal, and convolving.
But you never really convolve a signal by itself.
* I should never, ever, get on these boards before the caffeine has kicked in.  Yet, I regularly do.
** Well, probably quite intentionally similar characteristics, but -- caffeine.


Aliasing beating frequency As I understand aliasing is because you sample a high frequency with a low sample rate, and the sample points don't have enough information about the high frequency component. Here is an example of aliasing:

w = 5.0;(*signal frequency*)
dt = 0.66125;(*sampling rate*)
f[x_] := Sin[w x]
ls = Table[f[x], {x, dt, 200 dt, dt}];(*sample points*)

ListPlot[ls, PlotRange -> All, Joined -> True,
 DataRange -> {dt, 200 dt}]

we can see in the plot that there is beating.
And the Fourier transform gives only one frequency(except the symmetric negative frequency one).

DFT[A_, ht_] :=
 RotateRight[
  ht/Sqrt[2 \[Pi]]*
   Fourier[RotateLeft[A, Length[A]/2 - 1],
    FourierParameters -> {1, 1}], Length[A]/2 - 1]

ListPlot[Abs[DFT[ls, dt]]^2, PlotRange -> All, Joined -> True,
 DataRange -> {-1/dt/2, 1/dt/2}, Axes -> False, Frame -> True]

question:

Why there is the beating but only one peak in the Fourier transform? Since there is only one frequency in the sample, what is it beating against?
How to calculate the beating frequency and the aliasing frequency?

Update 1

As chirlu points out there is not beating, if it is the beating because of adjacent frequencies, we would expect to see two separated peaks if we increase the Fourier transform resolution. Here I increase the frequency resolution by including more sample points, and we can see it is a single peak, not adjacent peaks.

ls = Table[f[x], {x, dt, 800 dt, dt}];(*sample points*)ListPlot[ls,
 PlotRange -> All, Joined -> True, DataRange -> {dt, 200 dt}]

Update 2

Here is an example using sound:

Play[Sin[3000 t], {t, 0, 2}, SampleRate -> 6000]

Play[Sin[3000 t], {t, 0, 2}, SampleRate -> 950]

So use a lower sample rate, we can actually hear the beating!

Let’s look at an example. I’m going to use a sample rate of 50 throughout this post (i.e. the functions are sampled at 0, 0.02, 0.04, 0.06, …) and simple sine waves of the form $f_a(t) = \sin (2\pi\cdot a\cdot t)$, where $a$ is the frequency of the sine wave. For example, this is a plot of $f_3$:

The red impulses indicate the sample values. The green line, on the other hand, is the continuous function that is being sampled. The sample values all lie on this curve, and it goes through three complete cycles, as one would expect.

Aliasing

Because the sample rate is 50, anything above the Nyquist frequency of 25 is going to alias. As an example, here is a plot of $f_{53}$; mind the axes, I zoomed in to one tenth of the above, because otherwise the image would have become quite crowded.

There is slightly less than one sample per cycle, so the sampled values don’t follow the continuous waveform very closely. Actually, if all you saw were the sample values without the green curve, like so …

… you might be tempted to draw a much slower wave. Here is what happens when we plot not only $f_{53}$, but $f_3$ at the same time:

While the two sine waves are quite dissimilar in appearance, they happen to coincide at each and every sample point. There is absolutely no way to tell them apart given only the sample values. That’s aliasing. In fact, not only the samples of $f_{53}$ are exactly the same as those of $f_3$, but $f_{103}$, $f_{153}$, $f_{-47}$ etc. alias in the same way.

Put differently, due to sampling at a frequency of 50, the ranges 0–25, 50–75, 100–125, etc. collapse. What about the “upper half”, from the Nyquist frequency to the sampling frequency? It is mirrored; $f_{47}$ is the same as $-f_3$, and in general the frequency range 25–50 aliases to 25–0, with the samples inverted (or, equivalently, a phase shift by $\pi$).

This means, for your question, that the effect you are observing is not due to you using signals that exceed the Nyquist frequency; you can replace any higher-than-Nyquist-frequency signal with its alias in the Nyquist range and will get exactly the same results.

The “beating”

So what causes the apparent “beating”? Here is the very first plot of $f_3$ again, but this time it shows the sample values and linear interpolation between them (i.e. straight line segments from each sample point to the next), not the actual continuous waveform that was sampled:

The interpolated waveform is a bit edgy, but overall a good approximation of the original sine wave. Unfortunately, this is only because the frequency is very low relative to the sample rate. Here is $f_7$:

Urgh. It worsens as we get to higher frequencies, closer to the Nyquist frequency; this is $f_{23}$:

Unfortunately, even if we don’t draw those straight lines, our eye will try to connect the dots all on its own, and the impression is similar: It looks like a slow oscillation modulated onto a fast one. And that impression actually isn’t completely beside the point:

The same sample values, from $f_{23}$, can equally well explained by a ring modulation of $f_2$ and $f_{25}$ (with some phase adjustment), namely, $g(t) = -f_2(t)\cdot f_{25}(t+0.01)$. Thus, $g$ and $f_{23}$ are aliases, when sampled at a rate of 50. Here is a close-up again, showing that $g$ and $f_{23}$ in fact meet at 0, 0.02, 0.04 etc.:

A different way to write $g$ is $g(t) = 0.5 f_{23}(t) - 0.5 f_{27}(t)$, a sum of two sine waves of similar frequency, exactly the case where beating could be expected to happen. So it is beating after all, isn’t it? Not really, because $f_{27}$ is beyond the Nyquist frequency; as outlined above, it is an alias of $-f_{50-27}=-f_{23}$, and so $g(t)$ becomes $0.5 f_{23}(t) - 0.5 f_{27}(t) = 0.5 f_{23}(t) + 0.5 f_{23}(t) = f_{23}(t)$.

This is, within the baseband from 0 to 25, the only possible interpretation of those sample values.

Conclusion

Watch out for improper interpolation! When reconstructing a continuous signal from samples, the limited bandwidth must be taken into account (band-limited interpolation). Unfortunately, our eyes are not good at this; and many software waveform editors and resamplers aren’t, either.


Are the HOG features of the images good enough to do image classification alone? Why are their variances low? I've also asked this question here:https://cs.stackexchange.com/questions/77884/are-the-hog-features-of-the-images-good-enough-to-do-image-classification-alone

I've been working on classifying human poses like the pictures below. For that, I'm computing HOG features, and while in MATLAB, after using the command extractHOGfeatures, it's giving me 1764 HOG features, I also notice that maximum variance of all the 1764 HOG features are as low as 0.04 for a sample size of 2660, meaning the HOG features aren't probably changing much from sample-to-sample. My questions are:

1) why do you think we're getting low variances for HOG features? 2) what could we intuitively say about HOG features of these pose images? Are they good or interesting enough to classify images alone?

Definitely not sufficient to detect;

First HoG is not scale invariant - you change resize your image and handle this problem-.

Second HoG is not rotation invariant - you should include vary poses of action but this makes your solution more complex-

In my opinion; you should choose another kind of feature descriptors which can handle rotation and scale invariance. (SIFT, SURF, ORB ...)

Maybe you can use spatio-temporal information for this problem - It will be very useful-

Thanks.

And one more suggestion is about Neural Networks. If you have enough dataset for your problem, you can use "Neural Networks" which can be more roboust to other kind of feature based methods.


confusion in time invariance? I am trying to learn shift invariance(time invariance) and i came across a function/system

$$y(t)=3x(t)+2\cos(\pi t/3)$$

I am confused.If it was alone 3x,it should be definitely time invariant(i verified it through MATLAB)

But now it contains two terms, cos term in addition to 3x term,making function complex

Here is MATLAB code i tried to make for this system:

clc
clear all
close all
t=-5:.001:10;
x=heaviside(t)-heaviside(t-5);
y=3*x+2*cos(pi*t/3)
plot(t,y)

legend('y(t)')

figure
plot(t+3,y)

legend('y(t-3)')

figure
t=-5:.001:10;
x_shifted=heaviside(t-3)-heaviside(t-8);
y2=3*x_shifted+2*cos(pi*t/3);
plot(t,y2)

legend('S[x(t-3)]')

Its output shows that my system/function is not time invariant(please correct me if i am wrong)

The given system is not time invariant due to the term $\cos(\pi t / 3)$.

Given a system definition of the form:

$$y(t) = 3 x(t) + \cos( \pi t/3) $$ it's easy to see that

$$ y(t-d) = 3 x(t-d) + \cos(\pi (t-d) /3) \neq T\{x(t-d)\} = 3 x(t-d) + \cos(\pi t/3)$$

hence the system is time varying...