Orthornormal Wavelet scaling function coefficient calculation (filter design) I stumbled upon an equation system as shown in the image below while designing an orthonormal wavelet filter. I tried solving the equation (5.63,5.64,5.65 => 5.66) but no avail. I'm not sure if there's any property or parameterization method to solve this equation?

The equation below is taken from "Introduction to wavelets and wavelets transform A Primer" page 66

Add twice 5.65 to 5.64, and let $j = h(0)+h(2)$ and $k = h(1)+h(3)$; then we get

j + k = sqrt(2)
j*j + k*k = 1

which represents the intersection of a line with a circle, and we get

j = k = 1/sqrt(2)

Now subtract twice 5.65 from 5.64 and let $m = h(0)-h(2)$ and $n = h(1) - h(3)$; we get

m*m + n*n = 1

so for same angle $b$ say we get $m = cos(b)$ and $n = sin(b)$. We now have, for example

h(0) + h(2) = 1/sqrt(2)
h(0) - h(2) = cos(b)

which we can solve for $h(0)$ and $h(2)$. Finally shifting $b$ to get a gives the required expressions.


Meaningful output of FIR filter output I have implemented FIR filter using tapped delay line method. I start getting output as soon as first input sample is passed to it, I am wondering from which sample I will get the proper output(without transient outputs) from filter.

Does it depend on location of highest magnitude tap or it's just dependent on the length of filter.The number of filter coefficients I am using are N(odd) and symmetric.
for verifying I took the following filter taps as example :-

$F1$ $=$ $[1 0 0 0 0 0 0];$

$F2$ $=$ $[0 1 0 0 0 0 0];$

$F3$ $=$ $[0 0 1 0 0 0 0];$
and so on,if I pass sine wave using filter coefficients $F1$ then the first sample of the output is itself meaningful but if I use $F2$ and $F3$ then meaningful output is from second and third output samples respectively. I am also wondering how matlab command $filter()$ works.

If it's a linear phase FIR filter, then the inputs signals will be shifted (delayed) by an amount of group delay at the output.

For a linear phase FIR filter of length $L = 2K+1$ the group delay will be $N = K$ samples. For even length $L = 2K$ FIR filters it will be at $N = (L-1)/2$ ; half-sample position.

For non-linear phase FIR filters, group delay will be dependent on the specific frequency of the applied input signals.  It should be computed from

$$ \tau = - \frac{d \phi(\omega)}{d\omega} $$

where $\phi(\omega)$ is the phase response of the filter.

The matlab function filter(b,a,x) computes all the samples begining up to input signal length. So the first samples are transients.


Discrete (Inverse) Haar (2D) transform implementations I am interested in implementations of the Haar transform, and its inverse, such that when composed always result in the identity operation.  My inputs are discrete ordinals within a given range - and I need the output of the Haar transform to be similarly discrete ordinal (or fixed-point) numbers.

I've worked on the basis that if my input signal consists uniformly distributed samples on the range $0..(2^n-1)$ and, so does my transformed signal, then in principle, my original signal should be reproduced perfectly by the inverse haar.

I've dabbled with this (Matlab) implementation using GNU Octave:

http://people.sc.fsu.edu/~jburkardt/m_src/haar/haar.html

Which uses double precision values... forcing me to use round() to establish a discrete representation of the transformed data so as to fit the transformed representation into the same state-space as the original.  Unsurprisingly, I found I needed to use round() again to the output from the inverse transform.

It will probably come as no surprise that round(haar_2d_inverse(round(haar_2d(signal)))) is not quite the identity function for most values of signal... Anecdotally, there are typically a few out-by-one errors in the reconstructed signal - which seem roughly symmetric (+1 or -1 on a small number of the reconstructed samples) and this seems almost independent on the choice of n.

What I'd like to know is if there are 'better' implementations of haar_2d and haar_2d_inverse that work on ordinal samples?  Are the anomalies the consequence of the Haar implementation itself or the way I've applied round() at the intermediate stage?  If the latter, can I rectify this by scaling before rounding?

I think in 2D Haar, $\sqrt{2}$ is not required, if the transform is applied simultaneously in both directions. It'd be of the form of $\frac{a+b+c+d}{2}$, which requires only fixed point implementation with the Q point being chosen based on number of levels. I had used this to implement Haar in C for embedded systems, and it works without rounding/truncation errors.

Note:
My rep's not 50, so posting as answer. If anyone can move it to the comments it'd be good.


Measuring frequency response range of piezoelectric disc How can I measure the frequency response range of a piezoelectric disc below.

I don't have an oscilloscope but I do know how to program in matlab / octave I was thinking of creating a swept signal in matlab and playing it through my laptop or tablet and playing it through the disc.  But how do I measure the frequency response?

I am afraid that it is rather impossible without a proper hardware. Sweep sine is ok as a general method, but you would either need:

Reference transducer with known (preferably) linear frequency response. Then you can find the difference between those two.
Signal actuator with known frequency response, and then you can also find the difference between recorded signal and what it should be.

In any other case you will only measure the convolution of actuator response and response of piezoelectric element.

You should also understand that this disc is meant to be used as a contact device, not for recording of sound waves. There is a huge difference between those two signals recorded.

Lastly you can try to search for frequency characteristic of known model (presumably measured by someone). I did similar work back in time, but we used different type.


Is u[-n] defined from negative infinity to zero as it is folding of u[n] or from negative infinity to -1? Folding of unit step function and definition of its boundaries.

$u[n]$ is the unit-step function :

$$u[n] = \begin{cases} { 1 ~~~,~~~n \geq 0 \\ 0 ~~~,~~~ n < 0 } \end{cases}$$

then

$$u[-n] = \begin{cases} { 1 ~~~,~~~ -n \geq 0 \\ 0 ~~~,~~~ -n < 0 } \end{cases}$$

or stated as
$$u[-n] = \begin{cases} { 1 ~~~,~~~ n \leq 0 \\ 0 ~~~,~~~ n > 0 } \end{cases}$$

Hence, $u[-n]$ is nonzero for the range of $n \in (-\infty,0]$.


Number of taps needed in an FIR filter to remove DC An FIR filter is being used to remove the DC of an ECG. The sampling rate is 500Hz. How many taps would said FIR filter require (theoretically and practically) to filter out the DC.

P.S. the filter would be a high-pass filter and would not remove the 1Hz fundamental frequency component of the ECG.

See How many taps does an FIR filter need?

In your case you'd need more than 1000 taps depending on the allowable ripple, as your cut-off frequency is less than fs/500.

Alternatives :

use an IIR, a simple order-1 DC removal filter could work great
Average your signal and subtract the average in order to remove the DC
Rick Lyons proposes a clever implementation of an FIR DC removal filter using cascaded-integrator-comb filters here https://www.dsprelated.com/showarticle/58.php


What does the stability margin have to do with the quality of a filter? I designed a 2nd type Chebyshev bandpass filter with passband between 1 and 6 Hz and 30 dB attenuation in the stopband. When I looked at its Bode plot in Mathematica, I noticed "stability margins" (black lines) on the plot. I'm not familiar with this concept.

Does the stability margin have any impact on the stability of the IIR filter? For instance, will using the following filter on data (that is sampled at the appropriate sampling rate) cause stability issues?

The x-axis in both plots is the frequency (Hz) and the top plot is the magnitude response (dB) and the bottom is the unwrapped phase (radians).

Even though you might have designed your filter to be stable, Stability Margin is a concept which is used to estimate how far we are away from becoming Marginally stable/unstable both in terms of Gain and Phase.

This would mean then at Unity Gain(0 dB) any phase rotation of -180 degrees or ODD multiples of it is bound to make your system unstable and thus comes Phase Margin, the amount of rotation on the Pole Zero plot which would bring the Root Locus plot to intersect the point -1(or rather how far you're away from the 180 degree line on a Bode Plot when your Gain Magnitude is Unity).
Likewise we can state the same in terms of Gain. The amount of Gain that'd be necessary for the same Nyquist plot to expand and hit the point -1 (or how far you're away from the Unity magnitude gain when the phase response is -180 degrees).

Practical consequences of system implementation always require you to have some Margin because you're never guaranteed that the Gain and the Phase plots are exactly replicated when the system is actually working. So these parameters will tell you if you're safe or not. Theoretically, even being away from 1 degree away from -180 Degree line at Unity gain is stable; but practically there are things that can happen which will push it away easily.

You can learn a lot about these at

:http://en.wikipedia.org/wiki/Phase_margin
http://www.facstaff.bucknell.edu/mastascu/econtrolhtml/Freq/Nyquist4.html (a good simulation is shown here).

Stability Margin doesn't have any impact on the Stability of the Filter, it'll just indicate if your filter is Stable or not.

The Phase plot is actually shown in Degrees. Your filter looks like it is Marginally stable because at a few multiples of -180 degrees, the Gain is 0dB(oscillatory), and that too within the frequencies you've designed the filter to operate, it looks not promising at all. The Black Lines on the gain plot should appear at every multiples of -180 degrees. Also coming to the phase plot, there are no Black lines in this case because at least one Unity Gain crossover point has a Phase rotation that is an Integral odd multiple of -180 Degrees.


Are there any realtime voice activity detection (VAD) implementations available? I'm looking for a realtime VAD implementation like Discord/Mumble/Teamspeak have.

My goal is to have a system that let's me use the same VAD for all my VoIP applications, as well as something I can use to facilitate VAD for push to talk only applications.

Mumble is open source though I wasn't able to find what I needed when I had a look.

So, this is not an answer in the quality I'd like it to be, but the Mumble VAD is in AudioInput.c , as a quick browsing yielded:

speex_preprocess_ctl(sppPreprocess, SPEEX_PREPROCESS_GET_PROB, &prob);
fSpeechProb = static_cast<float>(prob) / 100.0f;

// clean microphone level: peak of filtered signal attenuated by AGC gain
dPeakCleanMic = qMax(dPeakSignal - gainValue, -96.0f);
float level = (g.s.vsVAD == Settings::SignalToNoise) ? fSpeechProb : (1.0f + dPeakCleanMic / 96.0f);

bool bIsSpeech = false;

if (level > g.s.fVADmax)
    bIsSpeech = true;

What's done here is that Speex, a FOSS voice codec and utility pack is used to get an estimate for the silence->voice transition. Whyever, Mumble converts that probability value (probabilities: $\in[0,1]$) to a percentage, and then compares it to the internal setting.

So the algorithm of interest here is definitely not within Mumble, but libspeexdsp!

Funnily, that's basically where it goes down in speexdsp's preprocess.dsp

/* FIXME: This VAD is a kludge */
st->speech_prob = Pframe;
if (st->vad_enabled)
{
    if (st->speech_prob > st->speech_prob_start || (st->was_speech && st->speech_prob > st->speech_prob_continue))
    {
        st->was_speech=1;
        return 1;
    } else
    {
        st->was_speech=0;
        return 0;
    }
} else {
    return 1;
}

So, this really just boils down to an SNR estimate of the current audio frame:

/* Speech probability of presence for the entire frame is based on the average filterbank a priori SNR */
Zframe = 0;
for (i=N;i<N+M;i++)
    Zframe = ADD32(Zframe, EXTEND32(st->zeta[i]));
Pframe = QCONST16(.1f,15)+MULT16_16_Q15(QCONST16(.899f,15),qcurve(DIV32_16(Zframe,st->nbands)));


What does "how to identify impulse response of a system?" mean? I don't know where to start from, I'm not familiar with math.
I'm trying to understand this question "how to identify impulse response of a system?".

As you can see as below figure, as I know, the input delta function and output is $h[n]$.

In this case, if I put impulse (function) signal in to above system, then I can get the impulse response. Then I think the output $h[n]$ and impulse response are the same, doesn't it?

Question 1: What is the meaning of "how to identify impulse response of a system"? Is this asking about convolution or auocorrelation? I can't understand what it does  mean.
Question 2: I want to know how to identify impulse response of a system?
Question 3: Why do we get the output of $y[n]$ when if we know that input $x[n]$, impulse response $h[n]$. I mean the output $y[n]$ is the same $h[n]$ isn't it? I'm very confused.

If you have any hint, please let me know.

Update

I found a question in my case as below,

  1.2.8 Consider a discrete-time LTI system described by the rule $$y(n)=x(n-5)+\frac 12 x(n-7).$$  What is the impulse response $h(n)$
  of the system ?

How to identify $h[n]$?

the answer :

$h[5]$ = 1+ 0 = 1

$h[7]$ = 0 + 1/2 = 1/2

. . .

Given a system, one often wants to be able to predict the output $y$, given an input $x$. This amounts to identifying a mathematical relation between all inputs and outputs, optimaly as $y = S(x)$. This can be a difficult task in general.

When the system is linear, and time-invariant (its properties do not change over time), then it is proven that it  suffices to know the system's response to a discrete unit pulse $\delta$, to be able to compute the output for any other input $x$, in the form of a convolution equation.

So suppose that your system outputs $h[n]$ when you input $\delta[n]$, then for any $x[n]$, the output will be:
$$y[n] = \sum_k x[k]h[n-k]\,.$$

To identify the impulse response of the system, you ought to provide the numbers, or a generic formula, that give the value for each $h[n]$.

You can try some exercices in Exercises in Signals, Systems, and Transforms, for  instance 1.2.4 and 1.2.7. You can also check the applet in the joy of convolution.

Since, in practice, it is impossible to generate a discrete pulse, there are other practical identification techniques for real-life systems, using random sequences or sine waves.


Why are Fourier analysis and transform only applicable for LTI systems? Why are Fourier analysis and transform only applicable for LTI systems?
What if the system is not LTI, won't Fourier analysis or transform be possible?

The family of Fourier Transforms are specificaly developed for analysing frequency contents of the signals for which there is no definition of linearity or time invariance. Hence we can define the Fourier transform of any signal, as long as it's integrable (i.e. stable).

On the other hand we can also define a Frequency Response $H(\omega)$ for a particular class of systems (those being LTI) when they are characterised by a unique signal denoted by $h(t)$ which is called as the impulse response of the LTI system.

The Fourier transform $H(\omega)$ of the impulse response $h(t)$ is important because the complex exponentials $e^{j\omega t}$ are eigenfunctions of such LTI systems and when the input is decomposed into a sum of such complex exponentials, via a Fourier transform, so will be the output, multiplied with a complex gain given by the Frequency Response $H(\omega)$ applied for each particular frequency. This gives a considerable ease for the development of the mathematical theory of signals and systems.

On the other hand consider the following nonlinear system defined as:
$$y(t) = x(t)^2$$
For which you can apply CTFT to both sides to derive the result in frequency:
$$Y(\omega) = \frac{1}{2\pi} X(\omega) \star X(\omega)$$

Also consider the following linear but Time-Varying system defined by:
$$y(t) = g(t)x(t)$$
For which you can also express the relation in Frequency domain via CTFT as:
$$Y(\omega) = \frac{1}{2\pi} G(\omega) \star X(\omega)$$

Which I hope convinces you about the scope of Fourier Transforms and their applications to systems other than LTI.


Recommended signal processing books? I am learning about signals in school but I found it very hard. Can someone help me by giving some references about any books?

We are studying Fourier transform, signal power, filters, digital and analog etc.

Without knowing what you find hard, it is difficult to recommend a specific book. Some are hard to read to. So, oldies but goldies, I suggest you to start with video lectures, and for instance by godfather A. V. Oppenheim

  Digital Signal Processing

and tell us where you have trouble (and why).


Trivial and non-trivial zeros I am new to DSP, and I'm self studying. Could someone please explain to me what do we mean by trivial and non-trivial zeros?

I wasn't familiar with that term in the context of signal processing. (Instead, I've seen the term being used in the context of the Riemann zeta function.) But I've found a document and this book where the term is used in a DSP context. The (obvious) definition is that trivial poles and zeros are the ones at the origin $z=0$ and at infinity $|z|=\infty$. They're called trivial because they don't affect the magnitude of the corresponding frequency response. Multiplying a given $\mathcal{Z}$-transform with $z$ (i.e. adding a pole at infinity and a zero at the origin) just advances the corresponding sequence by one sample, and multiplying by $z^{-1}$ (i.e., adding a pole at the origin and a zero at infinity) delays the corresponding sequence by one sample.

As an example, note that a causal FIR filter has as many poles as zeros, but all of them are at the origin $z=0$, i.e., they can't be used to shape the magnitude response. The same is true for an all-pole filter: all its zeros are at the origin, so they don't help to create a stop band (which would be one possible function of non-trivial zeros).

I think it's important to point out that the term trivial zero (or pole) is not a standard term used in DSP - as far as I know - but it appears to be an idiosyncratic use of a few authors.


Why does White Noise in images imply noise in adjacent pixels are independent? My professor said in the class that for "Additive White noise, the noise in pixels adjacent to each other are independent". How?

This is what I have got so far: White noise implies that PSD (Power Spectral Density) is flat which implies Covariance matrix is an impulse function i.e.

$$\begin{bmatrix} 1 & 0 & 0 & \cdots & 0 \\ 0 & 0 & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & 0\end{bmatrix}$$

How does this imply that noise at adjacent pixels is independent? At max, this says that noise is uncorrelated.

This question asks the same thing, though that part of the question is unanswered.

How did the professor state this? Did they say that the noise is iid?  That means that the samples are independent, identically distributed.  That means they are independent by definition (they are assumed to be so).

White noise, while a little more vague, also makes the assumption of independence. Again, this means that the noise is assumed to be independent.

However, if all that is said is that the samples are uncorrelated with the given covariance matrix, then you are correct: this does not imply independence.


Discrete Sinusoidal to State Space I'm looking to apply an optimal LQR filter to a discrete signal of the form

$x[n]=A\sin[\omega_0n + \phi]+ v[n]$

The amplitude $A$ and the phase $\phi$ are unknown variables I want to estimate using the filter, and $v[n]$ is an uncorrelated noise signal of variance ${\sigma_v}^2$.

I don't know how to build a state model to generate this sine wave, and proceed from there.

Thanks in advance for your help.

You may have a look on Paul Zarchan's lecture note about tracking a Sine Wave - Paul Zarchan - Fundamentals of Kalman Filtering: A Practical Approach - Tracking a Sine Wave.

You also may have a look on Tracking a Sine Wave with Kalman Filter - How to Account for Offset (DC Signal)?

Yet in practice, as @gsmafra wrote, using the DFT, since you know the frequency, is the best choice SNR wise.


How to filter this signal $x= \sin(2\pi 10 t)+\sin(2 \pi 40 t)+\sin(2 \pi 50 t)+\sin(2 \pi 60 t)$ using comb filters in time domain? Can anyone suggest a matlab implementation of this problem? I've done this uptill now,

fs = 4200;
t = 0.001:1/fs:1;
x= sin(2*pi*10*t)+sin(2*pi*40*t)+sin(2*pi*50*t)+sin(2*pi*60*t);
% figure,plot(t,x);
t1 = t - 0.05;
x1 = sin(2*pi*10*t1)+sin(2*pi*40*t1)+sin(2*pi*50*t1)+sin(2*pi*60*t1);
% figure,plot(t,x1);
yref = x - x1;
% plot(t1,yref);
N = 4096;
f = (0:N-1).*fs/N;
y_f = fft(yref,N);
plot(f,abs(y_f)),xlim([0 100]);
% freqz(y_f,1,1024,fs);

I am getting 10 and 50 Hz but from here how to remove 50 Hz?

I guess you are trying to design a nortch filter at 50Hz. You can try it in this way.

Generate a signal in frequency domain filter_f=ones(4200,1).
Then you set the points corresponding to 50Hz and -50Hz to 0.
Transform the signal to time domain.filter=ifftshift(ifft(filter_f)).
window the signal. filter_w=filter.*hamming(4200).
Pickthe center part points of filter_w. And that is the final nortch filter.


How to obtain impulse response from a plot of magnitude of frequency response? I have a plot of a magnitude response of a frequency response, how can I get the impulse response using only this plot and nothing else?

You state in your question that you have a "magnitude response". In that case, you cannot reconstruct the time-domain signal corresponding to the impulse response, because phase information is missing. In the plot, however, the vertical axis is indexed by $H(e^{j\omega})$ and not $|H(e^{j\omega})|$, so I suppose you meant "frequency response".

You can relate a system's frequency response to its respective impulse response via inverse DTFT:

$$H(e^{j\omega}) \xrightarrow{\mathscr{F^{-1}}} h(n)$$

Note that you have a typical ideal low-pass filter. That has a well-known DTFT inverse.


How to deal with modulator jitter? I have a question regarding an implementation of a modulator.
I guess some people here have an FPGA background.
I would like to implement a QPSK modulator of  31 Mbit bitrate using a 125 MHz clock.
I plan to generate the bitrate clock using MSB with a 32 bit phase accumulator also what can be called an NCO.
A 125 MHz clock has a 8 ns period.
A 31 MHz bit clock has a 32.26 ns period.
My bitrate clock will have an instantaneous period of 32 or 40 ns (but a mean of 32.26ns)
Is this a big issue to generate that kind of jitter on a modulator?
Is there a performance impact on a receiver?
Or maybe there is an another way implementing this?
Thanks

My bitrate clock will have an instantaneous period of 32 or 40 ns (but a mean of 32.26ns)

That won't work great. I'd hardly call a 20% too long period "jitter", anymore. You're simply not synthesizing a useful clock; you're average correct rate doesn't help a receiver much to get a lock on your symbol timing.
The "easy" solution is use the facilities of your FPGA to actually produce a 31 MHz clock.
If there's no clock synthesizer that can do that for you: You'll have to resample your signal internally. There's no need that your sample rate needs to be a multiple of your symbol rate, it's just easy to deal with.
But really, you can calculate what your transmit signal would look like at any rate (e.g. sampled at 62 MS/s) and then resample it to a rate that your system actually works at (e.g. 62.5 MS/s). There's no need for an actual 15.5 MHz clock to generate a baseband signal at 15.5 Msym/s QPSK to transport your 31 Mbit/s.
So, what you need is a rational resampler with a rate of 125/31 or 250/31 or so.
Your FPGA vendor quite likely has an IP core for that. If not, it's no magic, but will take some development time. You're looking for a polyphase implementation of a rational resampler if you need to implement that yourself, so that you can first decimate your signal by 31 and then interpolate by 125 (or 250, or so).


Doubts about Notch filter design I am following the process that is described in this question: Transfer function of second order notch filter , I want to create a notch filter with the band suppressed equal to $f_c = 4000$ Hz, so using $\omega_n= f_c / f_s $, ($f_s = 48000$), I obtained the $\omega_n = \frac{\pi}{6}$, then using the exact same formula, with $a =0.8$. The pole-zero graph I get has the zero in $1$ and a pole in $0.8$, is this correct?? I am getting the half of the filter since the filter is centered in $0$ and not in $4000$ Hz.

As far as I know it must be centered in the wn I have (based on $f =4000$ Hz), but I am not sure why it is centered in $0$, or how to center it in the desired frequency.
I get a pole zero graph like this one, with zeros in $1$ and poles in $0.8$.

Please advice

This Octave / Matlab code gives you a 2nd ord notch filter at $\omega_n = \pi/6$

r = 0.99;          % notch radius (closer to 1 stiffer)
wn = pi/6;         % notch radian frequency...

% Create the 2nd order NOTCH filter coefficients b() and a()
b = r*conv([1,  -exp(j*wn)],[1,  -exp(-j*wn)]);
a = conv([1,  -r*exp(j*wn)],[1,  -r*exp(-j*wn)]);

figure,freqz(b,a,2048);

with the following result:


Tips for programming scientific articles I have to reproduce simulations of scientific articles in Matlab.

Do you have any tips for programming scientific articles? What are your methods when you have to reproduce simulations of scientific articles?

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

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

Some of those are as follows:

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

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

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

This is very powerful for three reasons:

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

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

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

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

Hope this helps.


Difference between diversity combining and multiplexing How exactly are diversity combining and multiplexing different?

My current understanding is that, in certain types of diversity combining such as  maximal ratio combining, the same message is sent over different frequencies and then all the messages are collected and based on information of the SNR they are weighted and used as additional info to decode the message.

This seems to be the same concept involved in multiplexing. I am not quite understanding what is happening differently in multiplexing. Can someone please clarify?

diversity combining such as maximal ratio combining, the same message is sent over different frequencies

That is not inherently true.

Diversity combining makes use of different channels (not in the "TV channel" sense, but in the "everything that happens to a signal" sense).

Those channels might be created by over different frequencies all carrying the same signal, but often/usually, at the receiver, you use multiple e.g. antennas to observe the same signal.

The different signals are then weighted (often with a complex factor to reflect both phase and amplitude differences), and then added up, is what most commonly is referred to as diversity combining.

You could actually do the same with a single antenna and multiple identical signals mixed to different frequencies. However, that has little to no diversity gain: You actually need $N$ times the bandwidth of the original signal, so you could as well just transmit more data, add more redundancy etc.


Using MATLAB Function `mpiir_l2()`: Results Are Not IIR but FIR Filters - Why and How to Avoid This? I'm using mpiir_l2() from user Matt L's PhD thesis to design IIR filters. I set the number of numerator and denominator coefficients both to the same value (between $100$ and $300$). The maximum pole raduis is $0.98$. My observation is that the result is effectively a FIR filter, with the denominator coefficients allways being $[1, 0, 0, ...]$. Basically, the result is equal to that of the lslevin() function for least squares FIR design from the same source.

Clearly something in my filter spec makes the FIR solution have better errors than any IIR filter. The specification is the result of numerical optimization and looks pretty random to a humans eye, so I'm not sure what properties to look for.

Is there a way to get a true IIR filter from this function? I want to compare IIR anf FIR filters for the same specification and see what's better. Right now, I can't really compare when the results are basically the same.

EDIT: Here is a sample script and corresponding desired response. It seems like I have to use a third party hosting site for this, sorry.

close all;
clear all;

% loads the specification, weights and frequencies as three vectors
% the sepcification has a relevant band that is padded with zeros to both
% sides. There is a "don't care" gap between the relevant band and the
% zeros.
load 'desiredResponse.mat' % loads Ws, freqs, desiredResponse

% search area
allowedOrder = 1:1:3;       % numerator and denominator order
allowedShifts = 0:0.125:2;      % linear phase term

% other variables
Fs = 48000;
normalFreqs = freqs / Fs * 2 * pi; % normalized frequencies
poleRadius = 0.98;

% preallcoate for the search
lowestError = inf;
bestNum = [];
bestDenom = [];
errors = zeros(length(allowedOrder), length(allowedShifts));
orders = zeros(length(allowedOrder), length(allowedShifts));
shifts = zeros(length(allowedOrder), length(allowedShifts));

% search for best aproximation
for nl = 1:length(allowedOrder)
    ord = allowedOrder(nl);
    for ns = 1:length(allowedShifts)
        orders(nl, ns) = ord;
        shifts(nl, ns) = allowedShifts(ns);
        % add the linear phase shift
        phaseShiftTerm = exp(-1i * shifts(nl, ns) * (ord-1)/2/Fs * 2 * pi * freqs);
        D = desiredResponse .* phaseShiftTerm;

        % design the filter (change denominator order here, if you like)
        [num,denom] = mpiir_l2(ord, ord, normalFreqs, D, Ws, poleRadius);

        % calculate filter response
        H = freqz(num, denom, freqs, Fs);
        % calculate the weighted error
        errors(nl, ns) = sum(Ws .* abs(H - D).^2);

        if errors(nl, ns) < lowestError
            lowestError = errors(nl, ns);
            bestNum = num;
            bestDenom = denom;
        end
    end
end

I had a look at your specs, and I designed a few filters to see what is going on. First of all, we shouldn't expect that an IIR filter should perform much better than an FIR filter for the given specifications, because poles are mainly useful if sharp transitions from pass bands to stop bands must be realized, as is the case for frequency selective filters. In your case, we don't have sharp transitions, and adding poles is unlikely to improve the design substantially. Of course, theoretically, the approximation error should be decreased by adding poles (if the numerator order is not changed compared to the FIR solution), but the improvement over the FIR filter may be small. In practice, however, the resulting IIR filter could even be worse than the FIR filter with the same (numerator) order. There are two reasons for that: first, the design problem is much more difficult for the IIR filter and we might only find a local optimum, and, secondly, we need to restrict the pole radius to some value $r<1$ in order to get a stable filter that can be implemented. That restriction may increase the possible approximation error.

It indeed turns out that using an IIR filter does not improve things very much in your case. I used an FIR filter of length $N=100$ and an IIR filter with the same number of numerator coefficients, and $6$ additional poles away from the origin of the $z$-plane. The magnitudes of the frequency responses are shown in the top figure below (blue: desired, green: FIR, red: IIR). The IIR approximation is slightly better, but the difference might not be worth it. The bottom figure zooms in the lower frequencies, where the IIR filter approximates the desired response a bit better. The phase responses (not shown) are comparable for both filters. Note that for the FIR and the IIR designs I individually optimized the necessary linear phase term that was added to the given phase response.

The pole radii of the IIR filter are between $0.89$ and $0.98$ ($r=0.98$ was chosen as the maximum pole radius in the design process). So the designed IIR filter is really IIR, not just FIR with all its poles (approximately) at the origin. I do not know why you only obtain IIR designs that are actually FIR.


Different way to separate a particular frequency from a signal Let's say I have a signal which is composed of 2Hz,10Hz,17Hz,19Hz and 25Hz discrete time sinusoidal waves and I need to allow only 17Hz component to pass.

As far as I know, the standard way to approach the problem is implementing a Band Pass Filter.( Analog Prototype Filter -> Analog Filter -> Digital filter using Bilinear Transformation).

We can also implement DFT on the signal and use the information from the spectrum and then reconstruct the 17Hz component,right? If yes, which method is preferred and why?

If the frequencies to be rejected are the discrete frequencies given, and the sampling rate can be commensurate with these frequencies, then a moving average filter with the coefficients complex rotated to shift the response to 17 Hz can be done to reject all other integer frequencies as in the plot below.

This is the concept of "heterodyning the coefficients" rather than heterodyning the signal. When we heterodyne a signal, the signal is multiplied by a complex rotator ($e^{-j\omega t})$ which will in the process shift it in frequency according to the rate of rotation. In similar fashion, rotating your coefficients will shift the response of the filter. So you can either move your signal to the filter as in a downconversion followed by a low pass filter... or move your filter to the signal as in this case!
Here is the Matlab/Octave code that implemented this filter.
In this case the sampling rate is 26 Hz with a complex (I and Q) signal path.
z=exp(j*2*pi/26);
coeff=z.^([0:25]*17);
[h,w]=freqz(coeff,26,2048,'whole');
plot(w/pi*13,20*log10(abs(h)))
grid on
axis([0 26 -50 0])

(Note, this is the same frequency response as one bin in the DFT)
The coefficients themselves follow the pattern as shown in the following complex diagram, where the first 5 coefficients are labeled, and the rest can be followed along the trajectory lines shown. Observe that the magnitude of each coefficient is one (as in a moving average filter) but the phase is continuously rotated which causes the frequency response to shift at the rate of rotation.

Update:
As further suggested by @MarcusMüller in the comments, actually doing the signal homodyne followed by a low pass filter may likely be simpler in this case given that the OP has variable conditions for the signal of interest. Assuming the interference to be rejected is still discreet tones with an integer spacing relationship to the signal of interest (the signals need not be integer frequencies, but the multiplies of their spacing must be, as the nulls of a moving average filter will be located at n/T Hz where n is any integer and T is the time duration of the moving average in seconds).
So in this case we would first multiply the signal by $e^{-j2\pi 17 t}$ and then follow the (complex) output with a moving average filter (which means two filters, one for the I (in phase or real) path and one for the Q (quadrature or imaginary) path. In this case the moving average length as a minimum must be:
2Hz,10Hz,17Hz,19Hz and 25Hz
Move 17 to 0:
-12Hz, -7Hz, 0 Hz, 2Hz, 8Hz
And the least common multiple is 1 so a 1 second moving average is required with a commensurate sampling rate. Further as @MarcusMüller has suggested, this moving average structure could be implemented as a CIC decimating filter.  The output would be the "DC Equivalent" of the signal of interest, with the same amplitude and phase as the signal at 17 Hz. Further if the signal that is desired changes in frequency, the phase rotator can be changed in one place rather than updating the taps as I had done.
For further reading see:
Phase rotator: Numerically Controlled Oscillator (NCO) for phasor implementation?
CIC: CIC Cascaded Integrator-Comb spectrum
Note also to be considered as a filtering approach is the 2nd order resonator (or cascade of two of these structures), based on a rotated exponential averaging filter. (see FFT analysis for Vibration Signal). The image of this filter structure is duplicated below (see link for details on what the variables mean and how to compute them):


Why Is the Result of the Convolution of a Row and a Column a Two Dimensional Matrix? I'm trying to understand a book chapter on "Algorithms For Efficient Computation of Convolution" (Also on Doc Droid or Scribd) and I know that when calculating the convolution of an image represented by two dimensional matrix (array) with a mask matrix the result is the treated two dimensional image.

How to calculate the convolution of two signals represented by a row and a column?

The chapter I'm reading says that:

  Separable convolution kernel must fullfil the condition that its matrix has rank equal to one. In other words, all the rows must be linearly dependent. Why? Let us construct such a kernel. Given one row vector

  $$
\vec{u}=\left(u_1, u_2, u_3, \ldots, u_m\right)
$$
  and one column vector
  $$
\vec{v}^T=\left(v_1, v_2, v_3, \ldots, v_n\right)
$$
  let us convolve them together:
  $$
\vec{u}*\vec{v}=\left(u_1, u_2, u_3, \ldots, u_m\right)*\begin{pmatrix}v_1\\v_2\\v_3\\\vdots\\v_n\end{pmatrix}=\begin{pmatrix}u_1v_1&u_2v_1&u_3v_1&\cdots&u_mv_1\\
u_1v_2&u_2v_2&u_3v_2&\cdots&u_mv_2\\
u_1v_3&u_2v_3&u_3v_3&\cdots&u_mv_3\\
\vdots&\vdots&\vdots&\ddots&\vdots\\
u_1v_n&u_2v_n&u_3v_n&\cdots&u_mv_n\end{pmatrix}=A\tag{4}\\
$$

Following the wiki page on the subject, you are looking to create a 2 dimensional filter from two one-dimensional filters.

I found that the second example was more relevant for your question. I want to give a $\begin{bmatrix}1&2 &1 \end{bmatrix}$ weight column wise and on top of that give a $\begin{bmatrix}1\\2 \\1 \end{bmatrix}$ weight to each row. The result is the filter  $\begin{bmatrix}1&2 &1\\2&4&2 \\1&2 &1 \end{bmatrix}$ which take into consideration both your demands. By the way this spesific filter will smooth you image.

The simple fact that  this filter can be expressed as the outer product of two vectors renders this filter separable, Hence separable convolution. This post considers, also, non-separable filters.


The origin of constants in mel-scale formula To convert frequency to mel, we usually use formula:

$$\mathrm{mel}(f) = 1127 \ln \left(1  + \frac{f}{700}\right)$$

I wonder where are $1127$ and $700$ came from?

I've read the paper which Wikipedia think is the original of MFCC but I still didn't get the origin of those values. I hope somebody here can help me to understand the origin of those values.

The different mel scale formulas are supposed to approximate the human ear's critical bandwidths. Usually the formulas have this form:

$$m=C\log\left(1+\frac{f}{f_0}\right)\tag{1}$$

Below the frequency $f_0$ the mel scale changes approximately linearly with frequency, whereas above $f_0$ it changes logarithmically. This is the result of measurements. The difference between most formulas is the choice of the corner frequency $f_0$, which is usually chosen somewhere between $600\,\text{Hz}$ and $1000\,\text{Hz}$. So the number $700$ in your formula is the corner frequency $f_0$ where the scale changes from linear to logarithmic. The constant $C$ in $(1)$ is normally chosen such that $1000\,\text{Hz}$ correspond to $1000$ mel:

$$C=\frac{1000}{\log(1+1000/f_0)}\tag{2}$$

If we take the natural logarithm and with $f_0=700$ we get from $(2)$ $C=1127$. If instead we use the logarithm with base $10$, we obtain the other well-known constant $C=2595$.


minimizing the mean squared error? I'm reading this book "Applied Optimal Estimation" by A. Gelb. In the example 1.0-1, there are two measurements $z_{1}$ and $z_{2}$ with some noise for measuring $x$ as following
$$
z_{1} = x + v_{1} \\
z_{2} = x + v_{2}
$$
Now to design estimator for $x$, we need the estimated and the actual value of $x$, so we define the estimating error $\tilde{x}$ to be
$$
\tilde{x} = \hat{x} - x
$$
where $\hat{x} = k_{1}z_{1} + k_{2}z_{2}$. My question is that why $\hat{x}$ is the summation of the two measurements? Shouldn't it be one of the measurements or at least the summation divided by 2 so that we can make the subtraction between the ideal $x$ and observed $\hat{x}$ measurements.

Given measurements
$$\begin{align}
Z_1 &= x + N_1\\
Z_2 &= x + N_2
\end{align}$$
where $N_1$ and $N_2$ are independent zero-mean Gaussian random variables
with variances $\sigma_1^2$ and $\sigma_2^2$ respectively, it can be shown
that the minimum-mean-square-error (MMSE) estimate of $x$ in terms of $Z_1$ and
$Z_2$ is a linear function of $Z_1$ and $Z_2$.  Note that this linearity
is a consequence of the Gaussian assumption; for other kinds of distributions
that might be assumed for the noise variables $N_1$ and $N_2$, the MMSE
estimate is usually not a linear function.

So, let us consider the estimate $\hat{X} = k_0 + k_1Z_1 + k_2Z_2$ for the
unknown parameter $x$.  The mean-square error is thus
$$\begin{align}
E[(\hat{X}-x)^2] &=  E[(k_0 + k_1Z_1 + k_2Z_2 - x)^2]\\
&= E[(k_0 + (k_1+k_2-1)x + k_1N_1+k_2N_2)^2]\\
&= \operatorname{var}(k_1N_1+k_2N_2) + (k_0 + (k_1+k_2-1)x)^2
&{\scriptstyle{\text{because}~E[(Y+a)^2=\operatorname{var}(Y)+(\mu_Y+a)^2}}\\
&\geq \operatorname{var}(k_1N_1+k_2N_2)
\end{align}$$
where equality holds in that last inequality if we choose $k_0=0$ and $k_1+k_2=1$. With a slight change of notation to incorporate these
constraints into our calculations,
let us write
$$\begin{align}
E[(\hat{X}-x)^2] &= \operatorname{var}(\alpha N_1+(1-\alpha)N_2)\\
&= \alpha^2\sigma_1^2 + (1-\alpha)^2 \sigma_2^2
&{\scriptstyle{\text{because}~N_1, N_2 ~\text{are independent random variables}}}\\
&= \alpha^2(\sigma_1^2 + \sigma_2^2) -2\alpha\sigma_2^2 +  \sigma_2^2
\end{align}$$
where it is easily shown that the minimum value of that quadratic function
of $\alpha$ occurs when we choose
$$\alpha = \frac{\sigma_2^2}{\sigma_1^2 + \sigma_2^2} = k_1, \qquad
1-\alpha = \frac{\sigma_1^2}{\sigma_1^2 + \sigma_2^2} = k_2,$$
that is, when the two measurements are afflicted with Gaussian
noise of different variances (say with $\sigma_2^2 > \sigma_1^2$),
then the measurement with the smaller variance ($N_1$ in this
example) is given greater weight
$k_1=\frac{\sigma_2^2}{\sigma_1^2 + \sigma_2^2}$
than the measurement with the larger variance ($N_2$ in this case)
which has smaller weight $k_2 = \frac{\sigma_1^2}{\sigma_1^2 + \sigma_2^2}$.

When $\sigma_2^2 = \sigma_1^2$, the measurements are given equal weight
$\frac 12$ which is, according to the OP's insistence, the weight that
must be
assigned in all circumstances so as to minimize the mean-square error
because that is the only weighting that makes sense to the OP.


Amplitude problem by puting 5 FIR in a row I'm trying to use 5 FIR in a row to improve the S/N ratio on my signal, but I have some difficulties doing it.

In fact I have 2 signals, let's call them sig_I and sig_Q. I want to apply my 5 filter to them, and after that get the value PdB = 10*log( sqrt( sig_I² + sig_Q²) )

In order to do that the coefficients of the 4 first filters are given to me, and I have to be able to change the BW of the fifth one. It can be either 11 MHz or 60 MHz, but always with 32 coefficients.

At the end of this, I have an amplitude value, but my professor keeps telling me the amplitude is false, and that I should have about 7dB difference between the 2 cases (10*log(60/11) = 7.3 dB), but I can't find that.

A code example for that would be (on matlab) :

%Filter coefficients
b1 = [1 3 3 1]/8;
b2 = [1 3 3 1]/8;
b3 = [ 1 4 6 4 1]/16;
b4 = [ 1 5 10 10 5 1 ]/32;

fe = 2.4e9
bwFir = 5.5e6; % Can be 5.5e6 or 30e6
Wn = 2*bwFir/fe;
b5= fir1(32,Wn);

%convoluate signals with FIR
sig_I = conv(sig_I,b1);
sig_I = conv(sig_I,b2);
sig_I = conv(sig_I,b3);
sig_I = conv(sig_I,b4);
sig_I = conv(sig_I,b5);

sig_Q = conv(sig_Q,b1);
sig_Q = conv(sig_Q,b2);
sig_Q = conv(sig_Q,b3);
sig_Q = conv(sig_Q,b4);
sig_Q = conv(sig_Q,b5);

mod = 10.*log10(sqrt(I.*I + Q.*Q));

plot(mod);

But for me the value of "mod" is kinda the same with 11 or 60 MHz. Am I doing something wrong?

Look at the power spectral density of the signal prior to application of filter b5.

If you filtered white noise with two filters, one with a bandwidth of 11 MHz and the other with a bandwidth of 60 MHz, and both with the same passband gain, then indeed the difference between the power of the two output signals should be about 10*log10(60/11) or 7.3 dB.  But if, for example, your signal was already bandlimited to the 11 MHz bandwidth of filter 1, then there would be little difference between the output powers of the two filters.  A signal with a more complicated spectra (e.g., signal is within 11 MHz, but broadband noise exists) will be somewhere in between.

Also, if you are trying to compute the power of the signal in dB, you would not have a sqrt in your 10.*log10(sqrt(I.*I + Q.*Q)) equation.


Discretizing a Controller with the Backward Difference Method In the book "Modern Control Engineering" by Paraskevopoulos it's proved how to discretize a generic controller in the form:

$G(s)=\frac{Y(s)}{U(s)}=\frac{a}{s+a}$

where $a$ is a constant. Done this, the author says that it is possible to extend this procedure to a general case (without showing how), concluding that in order to discretize a controller $G(s)$ i have to simply substitute $s=\frac{1-z^{-1}}{T}$ where $T$ is the sampling period.

Since the author doesn't show how is possible to extend the procedure to a general controller. Can you explain to me ?

This is kind of hand-wavy, but you can look at this from two different perspectives:

One, you can look at $z^{-1}$ as a "back-step" operator; i.e. if $X(z) = \mathcal{Z}\lbrace x_n \rbrace$, then (with a few 'i's left undotted and 't's uncrossed) $\frac{X(z)}{z} = \mathcal{Z}\lbrace x_{n-1} \rbrace$.  You can also look at $s$ as a derivative operator: if $X(s) = \mathcal{Z}\lbrace x(t) \rbrace$, then $s X(s) = \mathcal{Z}\lbrace \frac{d}{dt}x(t) \rbrace$.

And -- the time derivative of $x(t)$ is approximately $\frac{d}{dt}x(t) \simeq \frac{x(t) - x(t - T)}{T}$.  If I throw caution to the wind and just sample the right side, that translates to, again roughly, $\frac{d}{dt}x(t) \simeq \frac{x_n - x_{n-1}}{T}$.  Now the Laplace of the left and the Z transform of the right gives $s X \simeq \frac{1 - z^{-1}}{T} X$.  Note that I have left off all pretense of a frequency domain variable on this -- I think that if there's an antithesis of "rigor" it must be "limpor", and we've just found it.

So you can take a system transfer function in $s$, expand it into a differential equation, then you can use the approximation $\frac{d}{dt}x(t) \simeq \frac{x(t) - x(t - T)}{T}$, then you can take the Z transform, then you turn the resulting difference equation into a transfer function in $z$ -- and you have just taken the very long way around to simply substituting $\frac{1 - z^{-1}}{T}$ for $s$.

The other way you can look at this is to observe that with a continuous-time system, the delay operator has a Laplace-domain transform: $\mathcal{L}\lbrace x(t - T)\rbrace = e^{-sT} X(s)$.  So you just set $z^{-1} = e^{-sT}$.  There's ways you can make this more rigorous, but compared to my above treatment this is already granite.

For $sT \ll 1$, $e^{-sT} \simeq 1 - sT$.  Substitute in $z^{-1} = e^{-sT}$ and solve for $s$ and you get $s \simeq \frac{1 - z^{-1}}{T}$.  Voila.  To get the first forward difference start with $z = e^{sT} \simeq 1 + sT$ and you find that $s \simeq \frac{z - 1}{T}$.  Use $z = e^{sT} = \frac{e^{\frac{sT}{2}}}{e^{-\frac{sT}{2}}}$ and you get the Tustin approxmation.

Note, when you actually go to do these, the word approximation.  You can generally get by with control systems because you need to sample way faster than the plant dynamics anyway, but -- approximation.


How to get an interpolation weight from a mathematical definition It was recently explained to me that a "Nearest neighbor" kernel for 1D interpolation can be implemented like this using NumPy

def nearest(delta):
    delta = abs(delta)
    if delta <= 0.5:
        return numpy.asarray([0,1])
    else:
return numpy.asarray([1,0])

Whereas the mathematical definition of nearest neighbor is

$h_{nn}(\delta) =
\begin{cases}
1 & \text{if}& -0.5 \le \delta < 0.5 \\
0 & otherwise \\
\end{cases}
$

Similarly linear interpolation, which in NumPy can be expressed as

def linear(delta):
    delta = abs(delta)
return [delta,1-delta]

But the mathematical definition for it goes

$h_{lin}(\delta) =
\begin{cases}
1-|\delta| &\text{if}&  0 \le |\delta| < 1 \\
0 & \text{if}& 1 \le |\delta|
\end{cases}
$

My question is how to form a kernel with weights out of these mathematical definitions, as the code belonging to them does not paint the same picture that the mathematical definitions paint.

I'll give a concrete example. Let's say you have a signal that was sampled every second. The frequency is $f = \frac{1}{1} = 1 Hz$.
Time: $T = \begin{bmatrix}0 & 1 & 2 & 3\end{bmatrix}$
Value: $X = \begin{bmatrix}1 & 2 & 3 & 4\end{bmatrix}$
We want to increase the frequency by $2$ i.e. $2 Hz$ (sample every 0.5 seconds).
Time: $\tilde{T} = \begin{bmatrix}0 & 0.5 & 1 & 1.5 & 2 & 2.5 & 3\end{bmatrix}$
Value: $\tilde{X} = \begin{bmatrix}1 & x_1 & 2 & x_2 & 3 & x_3 & 4\end{bmatrix}$
$x_1, x_2, x_3$ is determined by your interpolation functions $h_{nn}(\delta)$ and $h_{lin}(\delta)$. Here $\delta$ is time and both functions define intervals. Each value is given by $S(\delta) = \sum_{i=0}^{n-1} X_i \cdot h_{nn}(\delta-i)$. For nearest neighbor change the interval to $0 \leq \delta < 1$. Then $S(0) = S(0.5) = X_0 \cdot 1 + X_1 \cdot 0 + \cdots = X_0$. See also Cardinal B-splines.
While downsampling by $M$ requires strided convolution $$y[n] = \sum_k x[nM - k]h[k]$$
Upsampling needs fractionally strided convolution which is also called transposed convolution (see stackexchange)
$$y[j + nM] = \sum_k x[n-k]h[j+kM] \text{ and } j = 0, \dots, M-1$$
Transposed convolution with kernel size 3, stride 2 and padding 1 is equivalent to inserting 1 zero between inputs, pad by 1 and stride 1.
The kernel is $\begin{bmatrix}1 & 1 & 0\end{bmatrix}$ or $\begin{bmatrix}0 & 1 & 1\end{bmatrix}$ (either cross correlation or convolution) for nearest-neighbor interpolation (to double frequency):
from torch.nn import ConvTranspose1d
import torch
import numpy as np

def interpolate_nn(X):
  X = torch.from_numpy(X)
  with torch.no_grad():
    op = ConvTranspose1d(in_channels=1, out_channels=1,
                         kernel_size=3, stride=2,
                         bias=False, dilation=1, padding=1)
    op.weight.data = torch.tensor([0, 1, 1]).view(1, 1, -1).float()

    return op(X.view(1, 1, -1).float()).numpy().flatten()

X = np.array([1, 2, 3, 4])
print(interpolate_nn(X))

The result is [1. 1. 2. 2. 3. 3. 4.]
For linear interpolation use $\begin{bmatrix}0.5 & 1 & 0.5\end{bmatrix}$. The result is [1.  1.5 2.  2.5 3.  3.5 4. ]
Compare it with your $h_{lin}(\delta)$:
$\begin{align*}
S(0) &= X_0h_{lin}(0 - 0) + X_1h_{lin}(0 - 1) + \cdots = X_0(1 - |0|) = X_0\\
S(0.5) &= X_0h_{lin}(0.5 - 0) + X_1h_{lin}(0.5 - 1) + \cdots = 0.5X_0 + 0.5X_1\\
S(1) &= X_0h_{lin}(1 - 0) + X_1h_{lin}(1 - 1) + \cdots = 1X_1\\
\vdots
\end{align*}$


In more standard DSP terms, you have the following filter:

$$
y[n] = (1-a) x[n] + a y[n-1]
$$

where $x[n]$ and $y[n]$ are the input and output signals at time $n$ respectively.

The transfer function (which you didn't ask for) is:

$$
H(z) = \frac{1-a}{1 - az^{-1}}
$$

so here is your single pole, at $z=a$ in the complex plane. This filter is also known as exponential smoothing, exponential moving average (EMA), or exponentially weighted moving average (EWMA).

The infinite impulse response is $h[n] = (1-a) u[n] a^n$. In layman terms, when the input signal is 0 except for $x[0]=1$, the output signal is an exponential $(1-a) \times a^n$ starting at $n=0$.

What you want is the step response (i.e. what happens if the input signal is a constant $K$ starting at time $n=0$).

In this case, the output signal is the convolution $h \star Ku$ of the impulse response and the step signal. This is (for time $n \ge 0$):

$$
y[n] = \sum_{k=0}^{k=n} K \times h[k] = K(1-a) \sum_{k=0}^{k=n} a^k
= K(1-a)\frac{1-a^{n+1}}{1-a}
= K (1-a^{n+1})
$$

As the time $n$ grows, $a^{n+1}$ vanishes, and the step response grows monotonically to its limit value $y = K$, which is the value of the input signal.

This is what you get in your simulation.


Basic OFDM spectrum question As novice to OFDM, I have a question, related to modulator. Here, I present modulator block diagram, I found somewhere.

OFDM modulator AFAIK, final output signal has to have frequency spectrum consisted of $f_0, 2f_0, \ldots$ components positioned around $f_c$. Because inverse FFT produces complex numbers at it's output, my opinion is that signal at point A consists of pulses with different amplitudes having some frequency. If I am right so far, than the question is - what's the magic behind DAC block in order for him to have $f_0, 2f_0, \ldots$ components at his output B ?

The samples at the input to the IFFT block $X_0$, $X_1$, ... represent the magnitude and phase for each output frequency represented by the index. The inverse FFT translates this block to the equivalent time domain signal which has real and imaginary components (complex time domain representing your baseband signal). The dual DAC combined with a quadrature mixer converts this to a real signal centered at $f_c$ whose spectrum is identical to the baseband signal.

By declaring our data at the input to represent frequency with magnitude and phase represented by the constellation is central to the OFDM concept. You could imagine (properly) that each individual input $X_n$ to be a QPSK modulated waveform (such as your picture suggests), which is then up-converted onto its own sub-carrier, and banded together with the other sub-carriers with their own modulation, side by side in frequency, which is exactly what is occurring. The use of the IFFT just does this efficiently by putting each QPSK modulation onto it's appropriate sub-carrier by assigning it to an input of the IFFT (which must have frequency inputs), and allowing the IFFT to compute what the composite time domain signal of all of those sub-carriers would represent.


segmenting an EEG data I have an EEG labeled data, which is the data that used for training, And I want to segment those data based on the time of EEG signal (Time-based Epoching) as a preprocessing step,  based on the nature of EEG signals which is best:

To segment the data based on time domain, for example to take each 64 samples together.
Using the windowing functions like gaussian and hamming window. also if the best was this one,is it better to make the windows overlapped or not?

Overlapping segments don't cause any problem, unless your experimental design convinces you that such segments should not be allowed


wrong output using numpy fft on a tone I am getting the wrong output with the following simple implementation of Python fft.  Can someone please explain what I need to fix?

import numpy as np
import matplotlib.pyplot as plt

def tone(fs, frequency, dur=1, amp=50):
    '''
    test tone
    '''
    nsamples = dur * fs
    t = np.linspace(0, dur, nsamples, endpoint=False)

    tone = amp*np.sin(frequency*2*np.pi*t)
    return tone

fs = 16000
dur = 0.5
freq = 1500

stim = tone(fs, freq, dur=dur)
time = np.linspace(0, dur, num=fs*dur)

t = np.linspace(0, fs, num=fs*dur)
fdata = np.fft.fft(stim)

hwp = int(np.shape(fdata)[0]/2)
tplt = t[:hwp]
fplt = fdata[:hwp]

plt.figure(num=1)
plt.plot(time, stim, 'b-')
plt.show()

plt.figure(num=2)
plt.plot(tplt, np.abs(fplt.real), 'b-')
plt.xscale('log')
plt.show()

Removing

, endpoint = False

will fix the implementation though I do not understand why.
Perhaps someone can explain(?)

UPDATE: I suspected Spectral Leakage but that is not the case here. The reason for this is the OP is choosing to plot just the real component of the FFT, with a linear magnitude scale so the actual effect of the spectral leakage cannot be seen. (We actually get spectral leakage when we use endpoint = true!). If you change the plot from fplt.real to np.abs(fplt) you will see the true magnitude plot that must include the real and imaginary components!).

Here is the result with endpoint = true and plotting the dB of the magnitude.
plt.plot(tplt, 20*np.log10(np.abs(fplt)), 'b-'). The spectral leakage in this case is quite clear!

Here is the result with endpoint = false using
plt.plot(tplt, 20*np.log10(np.abs(fplt)), 'b-')

Here we see noise that is nearly 300 dB down (floating point precision).


Assume frequency response of x(t) is X(f). Find frequency response of x(-2t+4)? Frequency response of $x(t)$ is $X(f)$. Find frequency response of $x(-2t+4)?$
Here is my work using the fourier transform property:$
x(-2t+4) = x( -2(t-2) )$,
+ First, I apply the property,  $x(at)$  <---> $1/a \cdot X(f/a)$, by setting $a=2$, I have:
      $x(2t) <---> 1/2\cdot X(f/2)$
+ Second, I use the time-delay property,  $x(t-t0)$  <----> $X(f) \cdot e^{(-j \cdot 2\pi \cdot f \cdot t_0)}$, with $t_0=2$, I have:
      $x( 2(t-2) ) <---> 1/2 \cdot X(f/2) \cdot e^{(-j \cdot 2 \cdot \pi \cdot f \cdot 2)}$
+ Third, I use the rotation property, $x(-t)$   <-----> $X(-f)$, then I have,
     $x( -2(t-2) ) $<---> $1/2 \cdot X(-f/2)*e^{(-j \cdot 2 \cdot \pi \cdot f \cdot 2)}$
+ Lastly, because spectrum of FT is symmetric, means $X(-f)=X(f)$, so the final result is:
                        $1/2 \cdot X(f/2) \cdot e^{(-j \cdot 2 \cdot \pi \cdot f \cdot 2)}$
I am not sure what I did is wrong or not. Please help me point out and explain my wrong concepts or way?

The key is to apply 'shift and scale'. So for $x(-2t + 4)$, you would first do $x_1(t) = x(t+4)$, then $x_2(t) = x_1(-2t) = x(-2t + 4)$.

So $X_1(f) = e^{-j2\pi f(-4)}X(f)$ = $e^{j8\pi f}X(f)$.

Then $X_2(f) = \frac{1}{|-2|}X_1(f/-2)$ = $\frac{1}{2}e^{j-4\pi f}X(f/-2)$

In your answer how did you assume it is symmetric? Also, in your third step, for $x(-t)$, you should be getting $x(2(-t-2)$ instead of $x(-2(t-2))$. Your first 2 steps were correct but failed to use scaling factor -2 instead of 2 in the first step itself.


cceps (complex cepstrum) function in Matlab My doubt is in regards to the usage of the function cceps in Matlab.

The function is giving me the cepstrum for the entire speech signal of 3 or 4sec. This works well when I play back the signal after taking the inverse (using theicceps function). However, in order to analyze the given speech signal, we are supposed to only consider 20ms frame.

Is taking the complex cepstrum throughout the signal valid? Why is it working out here?

You should split the signal up into overlapping chunks of 20 to 30 ms. It is certainly valid to take the complex cepstrum on the 4 second signal. But, since the frequency characteristics of speech are constantly varying the complex cepstrum on the 4 second signal will average all of these together and you won't know what happened when. If you take it in each 20 ms window then you can see how the properties of the signal change over time.


Removing Noise from Dental Radiography I'm working on a project of applying Active Shape Model  to locate tooth in dental radiograph. For those familiar with the technique, I'm currently trying to sample along normal vectors for each landmark. The paper recommends to take derivatives of sampled pixels: "To reduce effects of global intensity changes we sample the derivative along the profile, rather than the absolute grey-level values."

So my problem is how to filter dental radiographs in the best way to prepare them for applying derivative operator. I'm currently using combination of median filter to remove most of what I think is quantum noise (mottle). It is followed by bilateral filter. Then I apply Scharr operator to compute the actual gradient which should be sampled.

The results are presented below:

The first image shows original data. In the second and third image, filtered data are presented, first as a magnitude of spectrum after FFT and then as a filtered image data. Fourth image show result of applying Scharr operator to third image.

My questions are:

Are there a well known approaches for reducing noise in dental radiograph that would differ from my approach?
What is causing the "smoky" appearance of the edges and "flat"(non-edge) areas? Is it some kind of leftover noise in the filtered image or is it inherent to the gradient operator? If it is indeed a noise, which filter would be most suitable to use? Median filter was good at removing small noisy blobs but large kernel causes the edges to blur out too much. So bilateral filter is used to filter out larger blobs and equalize color over area without harming edges, but it is not able to filter this smoky structure.
Is there a better option than Scharr operator to create gradient in this case?
Bonus: Would this be considered a good input for Active Shape Model? I'm not yet aware how robust they are.

As far as I understood, by image derivation you mean extracting edges. I would recommend to filter the image by a relatively large Gaussian filter. If computational cost of image derivation is uncritical to your work, I would recommend using canny edge detector. It is less sensitive to noise and does not fool by noise, and finds weak edges along with strong edges. Matlab instruction for that:

   [MinThresh MaxThresh]=[-0.3 0.5];
   EDGE_No_SMOKE=edge(im,'canny',[MinThresh MaxThresh]);

and the result is (I know it might not be the results you're seeking, however playing with Threshold variables and filter size will bring you desirable results):

Note that you do not see smoky effect anymore. Also about those wrong edges, you can remove them using image opening and closing techniques.


Decibels conversion into amplitude ratio I want to convert db into amplitude ratio.

I did some research but I only found online convertor, and I don't understand what's on Wikipedia very well...

Here is a little of backstory :
I am a programmer and I have to code an Equalizer; the problem is, I'm new to this DSP world.
I'm using a small documentation that an audio engineer in my company made a while ago and there are stuffs that I don't understand in it. (He's obviously not in the company anymore...)

In his documentation, he writes that the equalization is made step by step with an incrementation of ±0.1 dB.

So here are my questions :

What's the formula to convert dB into a percentage ?
On Wikipedia, there is 10 log10(x), but x is the power ratio, and I don't know how to get that either...
Does processing step by step is the best/only solution ?
Or can I calculate the ratio that I need each time (to limit the number of measures needed) ?

Thanks for you time.

You will want to work with what Wikipedia calls "field quantities."

Your engineer is telling you to work in steps of .1dB, which would be an amplitude ratio of .1=20log10(x).  To get x you do 10^(.1/20) = 1.01158

From x dB to ratio (amplitude):  10^(x/20)

From x ratio (amplitude) to dB:  20Log10(x)

If you can calculate the ratio, then you should be able to calculate the correction needed, round to the next .1dB step, and apply it in one step.


FFT for a specific frequency range. I would like to convert a signal to frequency domain. The desired frequency range is 0.1 Hz to 1 Hz and the frequency resolution is 0.01 Hz.

With the sampling rate of 30 Hz, FFT gives the frequency components up to 15 Hz.
Raising the sampling rate gives better frequency resolution. However, FFT gives wider frequency range. In my case, I just want 0.1 Hz to 1 Hz, FFT gives up to 15 Hz ( Extra computation ).

My question is, is there anyway standard way I can compute a frequency domain of a signal with specific frequency range and high resolution?

I think the best solution to your problem is to use the chirp-DFT. It's like a magnifying glass for a certain frequency range. It is more efficient than the direct implementation of the DFT (without FFT), because an FFT algorithm can be used with some appropriate pre- and post-processing. You basically need to modulate your signal with a chirp signal, then filter using an FFT, and then again chirp-modulate your signal to get the desired frequency response. See here and here for details on how to implement the chirp-DFT.


Estimating number of sinusoid from spectral estimation Suppose I am given a data X0 that have 50 realizations of a sum of an unknown number of sinusoids in additive white Gaussian noise.

I plotted the Spectral estimate as attached. How do I estimate the number of sinusoids in X0, their frequencies, and their relative amplitudes?
is it just the peak?

Edit:

n1 = 1; n2 = length(x)
Px = abs(fft(x(n1:n2),1024)).^2/(n2-n1+1);
Px(1)=Px(2);

The best approach and the details of how to do it depend greatly on the particulars of your situation.  If all your vectors are of the same length, and all your sinusoids have an integer number of cycles per vector (and therefore their frequencies correspond to a DFT bin) you can take the DFT of each vector.  Add the vectors coherently if the phases are the same in each random realization, add the DFT magnitudes if the phases are not the same.  Then set a threshold and find the points above the threshold.

If the frequencies are not on DFT bins but are spaced apart from each other, you need to modify the above approach to take account of the fact that the energy of each sinusoid is spread over multiple DFT bins.  You can still use threshold-based detection but will need to interpolate to find the frequencies.

If the frequencies could be spaced very closely together - which, from your plot, appears to be the case - the DFT-based approach is difficult.  The best approach may be a subspace-based technique like ESPRIT (Estimation of Signal Parameters Via Rotational Invariance Technique).  This involves estimating the covariance matrix and computing its eigendecomposition, then using the eigenvalues to estimate the power of each sinusoid and the eigenvectors to estimate the frequencies.  See, for examples: Richard Roy, Thomas Kailath: "Estimation of Signal Parameters via Rotational Invariance Techniques." IEEE Transactions on Acoustics, Speech, and Signal Processing.  Vol 37, No 7, July 1989.


Reason for bimodal behavior while low second fourier coefficient If I have a time series (for eg. for 23 timestamps) and if I plot it and see that it is bimodal, that means it might be having high value of second fourier coefficient (with frequency = 2). But when I do the fourier transform I find that the value of second fourier coefficient is very low.

Now, is the following reasoning meaningful for this?

It might be the case that although second coefficient has low value but other coefficients having high value are in phase with the second coefficient, hence adding to its amplitude, and that is the reason we are seeing bimodal behavior in original time series.

PS -
1) By second coefficient, I mean that when we do fourier transform we get a zeroth term for zero frequency, then first and then second.
2) By coefficient I mean the amplitude of the fourier coefficient.

For example, if I have a time series (bimodal i.e. with two peaks) and I do its fourier transform, 3rd term has an amplitude of x and phase of p1, where x is a very low value.

Now I also see that 2nd term has an amplitude of y and phase of p1 and 4th term has an amplitude of z and phase of p1. y and z have high values.

Now do you think that two peaks in the original time series can be attributed to the following reasoning -
Although 3rd term (for frequency=2, having two peaks) has low amplitude (x) but because 2nd and 4th terms have high amplitudes and same phase, they have added to the amplitude of 2nd term and hence we see two peaks in the original time series?

In general, instead of 2nd and 4th terms, it can be any coefficients having phase same as the 3rd term.

Two peaks does not imply a frequency of 2 - for example the peaks could be very narrow and/or their spacing could corresponding to a much smaller interval than half the time-series length.

There is no information about the overall shape of the time-series that can be drawn from the knowledge of the second DFT coefficient.

I have plotted here 15 time series (left) and the magnitude of their FFT (right). 5 unimodal, 5 bimodal, and 5 trimodal series.


Heterodyne modeling confusion in SDR I've been reading this lab sheet which explains the signal processing math of the RTL-SDR radio dongle. http://www.eas.uccs.edu/~mwickert/ece4670/lecture_notes/Lab6.pdf  In pages 5 and 6, the local oscillator (LO) is said to be modeled by a multiplication with the complex frequency $e^{-j2{\pi}f_{c}t}$ as shown in the image below taken from the lab sheet.
I understand how, mathematically, this multiplication with the negative frequency will shift the signal spectrum to the left by $f_c$ Hz, but what I don't understand is how a local oscillator can physically achieve this. The negative frequency in the behavioral model does not have a physical interpretation (does it?), and the physical LO signal in reality is just a sine wave with frequency $f_c$. Multiplying $sin(2{\pi}f_ct)$ with the radio signal of interest $s(t)$ will not remove the carrier frequency $f_c$ from the radio signal, for example if the radio signal is also $sin(2{\pi}f_ct)$, the resulting product will be a signal with frequency $2f_c$. Can someone explain what I am getting wrong?

A negative frequency is given by the complex signal $e^{-j2 \pi f_c t}$ for which the physical implementation requires two real signal paths, which is made clearer from Euler's formula:
$$e^{-j2 \pi f_c t} = \cos(2 \pi f_c t) + j \sin(2 \pi f_c t)$$
So if the Local Oscillator is generated with sine and cosine outputs, and the multiplier is a complex multiplier (in this case it may have a real input so would result in two real  multipliers) then the receiver can be implemented as described.
Note that a full complex multiplier requires four real multipliers:
$$(I_1+jQ_1)(I_2+jQ_2) = (I_1 I_2 - Q_1 Q_2) +  j(I_1Q_2- I_2Q_1)$$
In this case the input from the LNA is real, representing $I_1$ in the formula above, while the input from the Local Oscillator is complex, representing $I_2+jQ_2$ and the product is:
$$(I_1)(I_2+jQ_2) = (I_1 I_2) +  j(I_1Q_2)$$
With the complex output as indicated in the block diagram.
Note I detail this further and the implications of the real versus imaginary input after the LNA in this post: Frequency shifting of a quadrature mixed signal (for imaginary input the output of the LNA would be fed into a Hilbert Transform and then the multiplier would be a full complex multiplier using four real multipliers).


PSD estimate using FFT - scaling and units I know that similar questions have already been asked – I have read the answers, yet I am still not sure if I understand the topic properly at some points.

I want to use FFT to calculate PSD estimate. FFT gives me a set of complex numbers that I want to transform into the PSD values. I am aware of the fact that only the first $N/2 + 1$ values are useful, and that the other half are complex conjugates of the first half.

I know that I will need to take the magnitude squared of each (useful) complex number the FFT gave me in order to find the PSD values. I also found out that I should incorporate some scaling/normalizing of the result. It is at this point that I am unsure. In some of the replies to similar questions it is said that the magnitude squared should be multiplied by $1/N$ only, while in others it is said that it should be multiplied by $1 / \left(N\cdot F_s\right)$.  I have not been able to figure out which of these two should be used when. Can anyone explain?
The sampling frequency of my signal is $128\textrm{ Hz}$. The length of the signal to be transformed is $256$ samples. In this case I found out that the distance between the resultant values of the FFT (or PSD values) should be $128/256 = 0.5\textrm{ Hz}$. Is that correct?
The original signal is in $\mu V$. In what values will the result of the PSD be? Is it $\mu V^2$ or $\mu V^2 /\textrm{ Hz}$?
(There might be a connection here with what I am asking in 1), am I right?)

Have a look at this paper, section 9 goes over how to scale or normalize the FFT output. For the PSD you need to divide by the effective noise bandwidth (ENBW) which is

$$
ENBW = f_s\frac{S_2}{S_1^{2}}\\
S_1 = \sum_{j=0}^{N-1}w_j\\
S_2 = \sum_{j=0}^{N-1}w_j^{2}\\
where~w_j~is~the~jth~sample~of~the~window~function
$$


Is the Hilbert transform computable "losslessly"? The Hilbert transform is a procedure that can be used to perform a 90° phase shift on a signal.

Is the computation of (a discretized) Hilbert transform a "lossless" process? I.e. it merely phase shifts the signal, but does not alter it in any other way.

The reason for asking is that the following picture here looks like it introduces slight variation in the Hilbert transformed signal.

http://www.originlab.com/doc/Origin-Help/Hilbert-Transform

The loss you are seeing is due to the use of a finite length rectangular window.  If your signal is long enough that you can ignore the Hilbert filter transients at the both ends of the window (e.g. if it is infinite in length or so long that you don't care about the ends below some noise floor), then the process won't have the particular visible losses that your example points out.


What is the right way to extract particular frequency phase information from FFT? Signals acquisition:

Two proximity sensors are positioned Orthogonal to each other. These would be measuring the vibrations of a rotating shaft.

Problem:

Each of these 2 signals has two dominant frequency components (2441 Hz & 762.9 Hz ) in it. I don't know the phase and amplitude information. I'm interested in finding the phase of the signal at 2441 Hz and 762.9 Hz. I did the following process to achieve amplitudes and phases in MATLAB

Applying FFT to the raw signal and calculating Amplitude and Phase.

%fft for X signal
fftx=fft(X,NFFT);
% absolute value of the signal
Xn = abs(fftx);
A = mean(Xn);  % I took the mean to get single value.
phase = unwrap(angle(fftx));
phase_mean = mean(alpha);

Filtering raw signal to specific frequency prior to applying FFT using Butterworth filter. ( In this example I filtered to 2441 Hz)

[b,a] = butter(2, [ 2440/(fs/2),2441/(fs/2)],'bandpass');
X_filtered = filter(b,a,X);
% Applying FFT to filtered signal
fft_filtered = fft( X_filtered,NFFT);
Xn_filtered = abs(fft_filtered);
A_filtered = mean(Xn_filtered);
phase_filtered = unwrap(angle(fft_filtered));
phase_filtered_mean = mean(phase_filtered);

I did the same process as point 2 to filter to another frequency say 762.9 Hz.
My sampling frequency and other things are defined as follows.

L = length(X); % 50050 is length of signal
NFFT = 2^nextpow2(length(X)); % Zero padding to nearest N power 2
%Define frequency axis
fs = 1e7; % Sampling frequency
df = fs/NFFT; % frequency resolution
dt = 1/df; % time resolution
X_detrend = detrend(X,0); % Removing DC Offset

I've been reading a lot of other posts on the FFT phase and Amplitude.But I'm confused with most of them. In my case, I need to extract frequency particular information (amplitude and phase) to do further  calculation. I'm still unsure if using the filters is right way? Will FFT and filters alter the signal information in any way? Could someone enlighten in understanding the efficient way to extract amplitude and phase information from a signal of particular frequency component.

To estimate spectral component phase angle relative to shaft angle using an FFT of length N, first gather a sequence of N samples such that the shaft angle is known with accuracy at sample number N/2, then window with a Nuttal window and do an FFT shift before the FFT.  The  atan2() of the FFT bin result closest to the desired spectral component will produce an estimate of the phase relative to the center of the original sample window.  Adjust that phase by the known shaft angle at that N/2 sample point to get a spectral phase relative to shaft angle.

If you know the spectral frequency or use a frequency peak estimator, then the phase can also be interpolated between bins for increases phase estimation accuracy.  Interpolate phase by using a Sinc interpolation seperately on the real (even) and imaginary (odd) FFT result components before taking the arctangent.


Blurring transfer function of image I need help solving the following blurring function question.Assume an image $f(x,y)$ is moving in front of a camera so that $_0()$ and $_0()$ are the time-varying components of motion in the x- and y- directions respectively. The camera’s shutter opens at  = 0 and closes at  = . Assume that the shutter opening and closing operations are instantaneous. The blurred image (, ) captured by the camera is calculated as:

$$
 g(x,y) = \int_{0}^{T} f( x - x_0(t), y - y_0(t)) dt
$$

Calculate the blurring transfer function $H(u,v) = G(u,v) / F(u, v)$ in  frequency domain with respect to $_0()$ and $_0()$, where $(,)$ and $(,)$ are the Fourier transforms of $(,)$ and $(,)$ respectively and then calculate the blurring transfer function for the case where the image $(, )$ moves in only x- direction with a constant speed $\frac{a}{T}$.

Attempt: I really have no idea how to set up this blurring function.It is clear that the speed, $v$, is

$$
 v = x_0(t)i + y_0(t)j
$$

and also that

$$
 \mathcal F \{ f(x - x_0, y - y_0) \} = \text{exp}( (-2\pi i (x_0u + y_0v) ) )F(u,v)
$$

$$
G(u,v) = \mathcal F \{ g(x,y) \} = \mathcal F \left\{\int_{0}^{T} f( x - x_0(t), y - y_0(t))\right\}
$$

but I don't know how to proceed from there.Any help appreciated

Your start looks right, although I think your speed equation is wrong. To continue, the FT has no time dependance so you can take the integral out.

$$G(u, v) = \int_0^T \mathcal F\{f(x-x_0, y-y_0)\} dt $$

$$ = F(u,v) \int_0^T e^{-2\pi i(x_0u + y_0v)} dt$$

Divide by $F$ to get the transfer function.

If I understand the question correctly $x_0(t)$ and $y_0(t)$ are positions such that $x_0(t) = x_0(0) + v_x(t)t$ and similar for $y$. You can then put in an appropriate version of the speed given and do the integration.


What is the name of the following 1st order high-pass IIR filter? What is the name/family behind the following high-pass IIR filter:

$
T = 1/f_s \\
f_c = 5 \\
c_1 = \displaystyle \frac{1}{1 + \tan(\pi f_cT)}\\
c_2 = \displaystyle \frac{1 - \tan(\pi f_cT)}{1 + \tan(\pi f_cT)} \\
b = \left[c_1, -c_1\right] \\
a = \left[1, -c_2\right]
$

Where $f_s$ is the sampling rate and $f_c$ the cut-off frequency.

This is a first-order discrete-time Butterworth filter, which can be obtained via the bilinear transform from the following analog prototype filter (with normalized cut-off frequency):

$$H(s)=\frac{s}{1+s}\tag{1}$$

The bilinear transform replaces $s$ by $k\frac{z-1}{z+1}$, where $k$ is chosen such that the discrete-time filter has the specified cut-off frequency $f_c$, resulting in

$$k=\frac{1}{\tan(\pi f_c T)}\tag{2}$$

where $T=1/f_s$ is the sampling interval.

So if you replace $s$ in $(1)$ by

$$s=\frac{1}{\tan(\pi f_c T)}\frac{z-1}{z+1}\tag{3}$$

you obtain the corresponding discrete-time filter with the coefficients given in your question.

EDIT:

OK, here is how you obtain your final coefficients. Let $a=\tan(\pi f_c T)$. Plugging $(3)$ into $(1)$ then gives for the transfer function of the discrete-time filter

$$\begin{align}G(z)=\frac{\frac{1}{a}\frac{z-1}{z+1}}{1+\frac{1}{a}\frac{z-1}{z+1}}&=\frac{z-1}{a(z+1)+z-1}\\&=\frac{z-1}{(1+a)z-(1-a)}\\&=\frac{1}{1+a}\frac{z-1}{z-\frac{1-a}{1+a}}\end{align}\tag{4}$$

which is exactly in the given form with $c_1=1/(1+a)$ and $c_2=(1-a)/(1+a)$.


Is $\arctan(t)$ an energy or power signal? Is the function $\arctan(t)$  an energy or power signal?
I've tried working out the integrals for energy and power but I get stuck on both.

The typical inverse tangent function maps the input range of $ t \in (-\infty,\infty)$ into an output range of $(-\pi/2, \pi/2)$ as in the figure below:

Based on this, its values are bounded for all $t$. Yet, since the intergal of its square is unbounded, then it cannot be an energy signal;i.e.,

$$ \int_{-\infty}^{\infty} |\tan^{-1}(t)|^2 dt ~~~~~\to \infty $$

Moreover the average power of this signal is nonzero and finite, hence this is a power signal;i.e.,

$$ 0 < \lim_{T \to \infty} \frac{1}{T}\int_{-T/2}^{T/2} |\tan^{-1}(t)|^2 dt ~~~~~ < \infty $$

To evaluate the last integral, use the following MATLAB figure to see that

:

the square of the inverse tangent is always less than $(\pi/2)^2$ for all $t$. Hence to find an upper bound on the integral we can use $|\tan^{-1}(t)|^2 < (\pi/2)^2 = K$. Then the integral becomes

$$ 0 < \lim_{T \to \infty} \frac{1}{T}\int_{-T/2}^{T/2} |\tan^{-1}(t)|^2 dt < \lim_{T \to \infty} \frac{1}{T}\int_{-T/2}^{T/2} K dt = \lim_{T \to \infty} \frac{1}{T} K \cdot T = K ~~~~~ < \infty $$


QAM Constellation Construction - Basic 4/8 QAM I wanted to know how to construct the set of amplitudes/phase for QAM modulation. I saw that a QAM signal can be represented as:
$$
s_m(t)=A_{mi}g(t)\cos(2\pi f_c t)-A_{mq}g(t)\sin(2\pi f_c t)
$$
where $A_{mi}$ and $A_{mq}$are two amplitude levels and $g(t)$ is the pulse. The resultant amplitude is given by:
$$
A_{ra}=\sqrt{A_{mi}^2+A_{mq}^2}
$$

And the resultant phase is given by:
$$
\theta_m=\arctan\left(\frac{A_{mi}}{A_{mq}}\right)
$$

Given the above, what values should I take for the amplitude levels so that I will get the resultant amplitude and phase?

I tried generating the amplitude and phase values for simple 4 QAM by taking 2 different amplitude levels and 2 different phase levels like below:
$$A_{mi}=\left\{-1, 1\right\}, \quad A_{mq}=\left\{-1, 1\right\}$$

Then plotted the 4 possibilities (which seems pretty overlapping): Its clear that I need to take different amplitude levels for inphase and quadrature carriers. But is there any generalized formula for $A_{mi}$ and $A_{mq}$ so that I can generate my resultant amplitude and phase for QAM?

$$
\begin{array}{|c|c|c|c|}
\hline
A_{mi} & A_{mq} & A_{ra} & \theta_m\\\hline
-1 & -1 & \sqrt 2 & \frac \pi4\\\hline
-1 & 1 & \sqrt 2 & -\frac \pi4\\\hline
1 & -1 & \sqrt 2 & -\frac \pi4\\\hline
1 & 1 & \sqrt 2 & \frac \pi4\\\hline
\end{array}
$$

You can choose any mapping that you desire from your unique data sequences to the symbol locations on the QAM constellation. It is common and usually best to Gray-code the mapping such that only one bit changes between adjacent constellation locations (those with minimum Euclidean distance). This way if there is a symbol error due to noise crossing an immediate decision threshold, there will only be one bit error.

See figure below showing the common mapping for QAM using Gray Coding from Wikipedia under "Gray Coding". Notice that the adjacent constellation points only differ by one bit. Under AWGN conditions at threshold error conditions with equi-probable symbols only one bit error would result with this mapping approach (the optimum error boundaries with equi-probable constellation points are mid way between the points shown).

For Gray-coded QPSK the mapping is quite trivial as illustrated by the figure from this link:

And here is a link to a simulation showing the significance in BER performance when Gray-coding QPSK.

The amplitude and phase for each possible symbol is directly taken from the constellation's I and Q values as shown in the diagrams above, and as traditionally done for any complex I, Q pair:

The amplitude is

$$A= \sqrt{I^2+Q^2}$$

And the phase is:

$$\phi = \mathrm{atan2}(Q,I)$$


Naive Kalman filter for 3D position I looked at posts that discusses 3D kalman filter.
Kalman Filter to estimate 3D position of a node
Help with Kalman Filter implementation for estimating 3D position

Both from my understanding, both person were attempting the same thing, which is calculating 3D position based on accelerator readings. Yet their process for setting up the Kalman filter is completely different.

My own method is also different

These are my equation based on what I learned here, it treats acceleration as the U vector http://www.cl.cam.ac.uk/~rmf25/papers/Understanding%20the%20Basis%20of%20the%20Kalman%20Filter.pdf

Xk = AXk-1 + BUt
x_curr = x_prev + x_vel_prev * dt + 0.5*dt^2*x_acc
x_vel_curr = x_vel_prev + dt*x_acc

The sescond post seemed to have dropped  the U vector completely similiar to this example http://campar.in.tum.de/Chair/KalmanFilter#Literature.

While the first post has a different process equation based on wiki and it uses Ga instead of Bu https://en.wikipedia.org/wiki/Kalman_filter#Example_application.2C_technical:

X = FX_prev + Ga

This is what I did and it make sense I think:

X = [x,y,z,dx,dy,dz]

U = [ddx, ddy, ddz]

A = np.matrix([[1, 0, 0, dt, 0,  0],
               [0, 1, 0, 0,  dt, 0],
               [0, 0, 1, 0,  0,  dt],
               [0, 0, 0, 1,  0,  0],
               [0, 0, 0, 0,  1,  0],
               [0, 0, 0, 0,  0,  1]],
             )
B = np.matrix([[0.5*(dt**2), 0,           0],
               [0,           0.5*(dt**2), 0],
               [0,           0,           0.5*(dt**2)],
               [dt,          0,           0],
               [0,           dt,          0],
               [0,           0,           dt]])

H = np.matrix([[1,0,0,0,0,0],
              [0,1,0,0,0,0],
              [0,0,1,0,0,0]])
P = np.matrix([[10,0,0,0,0,0],
              [0,10,0,0,0,0],
              [0,0,10,0,0,0],
              [0,0,0,100,0,0],
              [0,0,0,0,100,0],
              [0,0,0,0,0,100]])

Q = np.matrix([[1,0,0,0,0,0],
              [0,1,0,0,0,0],
              [0,0,1,0,0,0],
              [0,0,0,1,0,0],
              [0,0,0,0,1,0],
              [0,0,0,0,0,1]])

R = np.matrix([[1,0,0],
              [0,1,0],
              [0,0,1]])

What I dont understand is, how come the previous 2 posts have completely different matrix? Is there a reason why? And is my approach correct?

Any input and suggestion is greatly appreciated

The thing you should keep in the back of your mind is that R Kalman published his seminal paper in 1962 (or there about).  He was a Control Engineer, interested in real systems at a time prior to widely available high fidelity numerical multi-physics modeling. The kinematic models you contrast, are derived from
$$
F=ma =m \frac{dv}{dt}.
$$
and $m$ is often a constant point mass.
  Planes and rockets consume mass for propulsion so we can see early on that there will be neglected terms.  $F$ should also include aerodynamic drag and for laminar flow, drag is roughly proportional $v^2$ which is not even close to the Navier-Stokes  equations.  So even a simple drag model is nonlinear and for a real time online algorithm a linear state space model is all he really had.  The solution to these oversimplifications is noise and small step sizes and prediction feedback.

Noise can also be actual random phenomena.  So a lot of non trivial modeling errors are lumped together in
$$
\frac{d}{dt} \mathbf{x}(t)= \mathbf{A}\mathbf{x}(t) + \mathbf{B} \mathbf{u}(t)
$$
where $\mathbf{u}(t)$ can be control forces, random forces, or psuedo modeling errors that we treat as noise. We also have a measurement set of equations that can also include noise, real or conceptual.
You can stick with this form and use a Kalman-Bucy Filter or you can convert the continuous state space model to a difference equation which is typically used in what is a Kalman Filter. To do that you need to integrate the equation just given.  I'm not going to go in detail because there are books and websites that show this but one start with calculating the matrix exponnential
$$\mathbf{F}= \mathbf{e}^{\mathbf{A} \Delta t} $$.  Some of your examples calculate the matrix exponential and some don't. You can often obtain a symbolic solutution for the matrix exponential or you can solve it numerically.

Getting back to your original question is that these models are simplifications of simplifications.  There may be poor choices taken but there really isn't a single correct discrete time 3D state space model. some may have $x_{k+1}=x_{k} + \Delta t v_{k}$ and some might include $\frac{(\Delta t)^2}{2} a_k$.  Real objects aren't point masses.  Angular orientation matters.  A lot of what is missing is treated as noise over small time steps.  When the nonlinearity is more severe, we switch to the extended Kalman Filter.  In the future, the state models will be real time multi-physics numerical routines.

The best model is the one that works best with you data.


down sampling and a factor of 2 without loss of information I am trying to understand the limits of down sampling on a signal, i.e. how far can I go without loss of information or overlapping. It says in the solution manual of the Oppenhiem and Schafer book (on a question where down sampling is 2) that

"There is no loss of information if X(e^(jw/2) ) and X( e^(j(w/2) - pi) ) do not overlap."

Does this mean that my original signal is periodic and when I down sample means that my frequency response X(e^jw) becomes X(e^(jw/2) ) ?

Please help,

If you down-sample without filtering, you will fold spectrum above half the new sample rate together with spectrum below half the new sample rate.  If both those spectral bands (above and below) are sparse enough (big gaps in the spectral content) in the right places, then it is possible that when folded, there will be no overlap of the folded non-zero high band spectrum with the non-zero low band spectrum, thus no loss of information due to aliasing.  If you know where the spectral gaps were, you could unfold the data (with some effort).

In the more typical case where those two spectrum bands (upper and lower) are not sparse enough and thus do overlap, the two spectrums be mixed (like mixing paint) when downsampling, thus incur an informationally lossy process.


Convolution by differentiation property of convolution This is the question

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

Actual answer should come like this

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

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

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

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

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

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

the left side evaluates as;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Sorry for the confusion...


A reference request for zero-order hold (ZOH) From Wikipedia:

  A zero-order hold reconstructs the following continuous-time waveform
  from a sample sequence $x[n]$, assuming one sample per time interval
  $T$: $$x_{\mathrm{ZOH}}(t)\,= \sum_{n=-\infty}^{\infty} x[n]\cdot
 \mathrm{rect} \left(\frac{t-T/2 -nT}{T} \right) \ $$ where
  $\mathrm{rect}()$  is the rectangular function.

I'd like to find a good bibliographic reference for it.

Thanks in advance.

I'd personally start with "Signals and Systems" by Oppenheim and Willsky. I have the second edition, and it has a good description of zero-order hold. This book is a must-have if you're interested in signals and systems. There is also a version that focuses on discrete signals. In addition, this book is probably available at most engineering libraries.