M-Array DPSK decoding and encoding So i'm trying to figure how how to do the encoding and decoding of an M-array Differential Phase Shift Keying (DPSK) signal. I understand how to do it for binary DPSK, but on higher ordered signals I'm missing something. Here's how I would do it for binary DPSK (M=2).

The steps should be as follows: get Data sequence -> Differential encoding -> assign phase for Tx -> (on the receiver side) do Differential decoding - assign bits from phase. The following assumes an initial bit of 1, where 1 corresponds to a phase of $\pi$

$\;\;\;\;$Data Sequence $b$ =   | 1 0 1 1 0 0 1 0 1 0

$\;\;\;\;$Diff Encoding $a$ = 1 | 0 0 1 0 0 0 1 1 0 0  -> this opperation should be $a_k = a_{k-1}\oplus b_k$

$\;\;\;\;$Assign phase $\phi$= $\pi$ | 0 0 $\pi$ 0 0 0 $\pi$ $\pi$ 0 0

$\;\;\;\;$Diff Decoding $\rho$= $\ $$\ $$\ $| $\pi$ 0 $\pi$ $\pi$ 0 0 $\pi$ 0 $\pi$ 0 -> where $\rho = \phi_n - \phi_{n-1}$

$\;\;\;\;$Phase to bit = $\;\;\;\;\;\;\;$| 1 0 1 1 0 0 1 0 1 0

So I end up with my original sequence. but if I try to do this with M=4, I don't end up with the same sequence. Lets assume for M = 4 that 00 = 0, 01 = $\frac{\pi}{2}$, 11 = $\pi$ , and 10 = $\frac{3\pi}{2}$, and that we start with initial data 11

$\;\;\;\;$Data Sequence $b$ = $\ $  | 10 11 00 10 10

$\;\;\;\;$Diff Encoding $a$ = 11 | 00 11 11 01 11  -> this opperation should be $a_k = a_{k-1}\oplus b_k$ = ($a_{k-1} + b_k$) mod M

$\;\;\;\;$Assign phase $\phi$= $\ $ $\pi$ | 0 $\pi$ $\pi$ $\frac{\pi}{2}$ $\pi$

$\;\;\;\;$Diff Decoding $\rho$=$\ $$\ $$\ $$\ $$\ $ | $\pi$ $\pi$ 0 $\frac{3\pi}{2}$ $\frac{\pi}{2}$ ->  where $\rho = \phi_n - \phi_{n-1}$

$\;\;\;\;$Phase to bit =$\, $ $\, $$\, $$\, $$\, $$\, $$\, $$\, $$\, $$\, $ | 11 11 00 10 01

So here I don't end up with the same sequence. I'm not quite sure where I'm messing up there. If feel like if I can do this for M=4 I can do this for any size M. Can someone help me see where I'm going wrong? Thanks for your time!

As vaz's answer suggests, your differential encoding is incorrect. In QPSK
regarded as a phase-shift-keyed signaling scheme as opposed to
a 4-QAM scheme (see this previous answer of mine for some details), the input dibit $(b_I,b_Q)$ is treated as being the Gray-coded integer representation of
the integers $\{0, 1, 2, 3\}$ and the phase of the RF carrier in the
corresponding signaling interval is $0$ or $\pi/2$ or $\pi$ or $3\pi/2$. In tabular form, we have
$$
\begin{array}{|c|c|c|c|c|c|}
\hline
(b_I,b_Q) & \text{normal value} ~k & \text{Gray code value} ~\ell  &\text{phase-modulated signal}\\
\hline
(0,0) & 0 & 0
& \sqrt{2}\cos\left(2\pi f_c t - 0\frac{\pi}{2}\right)\\
(0,1) & 1 & 1
& \sqrt{2}\cos\left(2\pi f_c t - 1\frac{\pi}{2}\right)\\
(1,1) & 3 & 2
& \sqrt{2}\cos\left(2\pi f_c t - 2\frac{\pi}{2}\right)\\
(1,0) & 2 & 3
& \sqrt{2}\cos\left(2\pi f_c t - 3\frac{\pi}{2}\right)\\
\hline
\end{array}
$$
That is, we can say that the QPSK modulator input
$(b_I,b_Q)$ is regarded as the Gray code representation
of the integer $\ell \in \{0,1,2,3\}$ and the modulator produces the
output
$$\sqrt{2}\cos\left(2\pi f_c t - \ell\frac{\pi}{2}\right).$$

In differential QPSK (DQPSK), the input dibit Gray-code value
is subtracted (differential, get it?) from the previous transmitted dibit Gray-code value (left-handed folks can use addition if they prefer) to get
the next transmitted dibit Gray code value.  That is, with $\ell$
continuing to stand for Gray code value, we have that
$$a_{\ell}[n] \equiv a_{\ell}[n-1] \mp b_{\ell}[n] \bmod 4. \tag{1}$$
At the demodulator, we have
$$\hat{b}_{\ell}[n]\equiv \hat{a}_{\ell}[n-1] \mp \hat{a}_{\ell}[n] \bmod 4 \tag{2}$$ where the carets denote the demodulated symbols and
it is the job of the project manager to make sure that
the modulator designer and the demodulator designer both use the same
sign ($+$ or $-$) in $(1)$ and $(2)$.  In well-designed DQPSK systems,
there is no explicit (two-bit binary) arithmetic used in $(2)$,
and the DQPSK demodulator determines $b_{\ell}[n]$ from soft-decision
data or more usually, the raw matched-filter/correlator outputs in
the I and Q branches at two successive sampling instants.

Finally, you can apply the same idea to $M$-PSK with $M = 2^n > 4$
with the successive phase shifts being labeled with the Gray-coded representations of the integers $\{0, 1, 2, \ldots, 2^n-1\}$.
Equations $(1)$ and $(2)$ above get modified to be arithmetic modulo
$M$ but the same ideas work. The receivers are more complicated, though.


What is the correct way to apply multiple IIR filters to a signal in MATLAB I am attempting to apply multiple IIR filters (one band-pass, one high-pass, and one low-pass;each with Fs of 44.1 kHz ) to a signal in MATLAB by using the filter function.  When I apply one signal and sound it, I can hear it. But as I start to apply the others, I can no longer hear the signal. I have the following code:

 audio = handles.filename;
 [x,Fs] = audioread(audio);
 lowpass = IIR_Lowpass_Butterworth_CFL1000;
 highpass = IIR_Highpass_Butterworth_CFH10000;
 bandpass = IIR_Bandpass_Butterworth_CFL1000_CFH10000;

 y = filter(lowpass,x);
 z = filter(highpass,y);
 q = filter(bandpass,z);

 sound(q);

Am I applying the filters incorrectly?

If you filter out all frequencies above a certain cut-off frequency by a low pass filter, and if the filter is good, then there's hardly anything left above that cut-off frequency. If you then filter out everything below a certain cut-off frequency using a high pass filter, and if that cut-off frequency is greater than the cut-off of the low pass, then you eliminate the rest. There's not much left to do for the band pass filter afterwards, because you've basically suppressed your whole signal. That's why you don't hear anything.

What you probably wanted to do was apply the filters separately, i.e. create one low pass output, one high pass output, and one band pass output. In this case you have to filter the original signal x with the three different filters to generate three different output signals.

yl = filter(lowpass,x);
yh = filter(highpass,x);
yb = filter(bandpass,x);


Vertical lines noise in spectrogram I have started looking into DSP and made an implementation of short-time fourier transform, using kissFFT to perform the FFT. I then render a spectrogram.

An example image that I have generated can be seen below. On top is the spectrogram of a music file from my implementation and under it is the spectrogram of the same section generated from audacity.

My concern is the vertical lines that show up from my implementation, which are most visible in the higher frequencies. They may interfere with other stuff I will try to do later on. Has anyone ever had the same problem or know the cause? How can I get rid of them?

Additional information:
I use hamming window with 50% overlap (using different windows has not helped)

After wasting forever on this problem, I found that it was just a question of window function. Doh!

Window functions that are not zero-ended produce these vertical lines. For example, Rectangular Window, Hamming Window, Gaussian Window (with low sigma) produce these lines, and Barlett Window, Hanning Window, Blackmann Window, Welch Window, don't.

I had tried other window functions but maybe during my tests, I made a coding mistake so I didn't find the fix. Moral of the story : do tests carefully!


How to read video in octave? I am working on a color transformation project using video.Can any one please tell me how to read video in octave.I have see aviread can be used to video frame by frame but i want to read the entire video file.

There is a package from octave. I am using debian jessie and this version of octave-video producing errors and can not be installed. Also, there is no debian package with name octave-video.


How to approximate the sample rate? I have a simple program that captures audio from an audio device. I have configured a nominal sample rate of 48000Hz and a buffer size of 1 millisecond. The audio sub system should execute my capture code with approximately 1000Hz.

My assumption is that the audio clock (LRCLK) is not running at exactly 48000Hz. Using this 1kHz-callback from the audio system I want to approximate the actual sampling frequency relativ to another clock. Currently I use the system clock of my workstation as a reference. I plan to use a network clock.

For my very first experiment I used a simple moving average with a period of 1000 but there is still very much jitter. Sometimes the average jumps by 2Hz which seems much to me.

What is the proper way to go here?

Many thanks!

In my experience, sound card sampling clocks vary by as much as 10 Hz from one card to another, even between computers of the same make and model. I have never seen much jitter, but I haven't especifically looked for it.

Note, however, that it is essentially impossible to know the true sampling rate, unless you have access to an atomic clock. What you can do is approximate the difference relative to a clock that you designate as your reference.

A simple experiment is to generate a single-tone signal, then receive it and find its frequency. Say your reference is the transmitter and it's operating at 48,000 samples per second. Feed the transmitter a sequence 0, 1, 0, -1, 0, .... This will produce a sine wave of nominal frequency 12,000 Hz. In the receiver, you sample this signal for a relatively long time (to have good frequency resolution) and use the DFT to determine its frequency to be, say, $f_2$. Then, you can calculate $$\Delta=\frac{12000}{f_2}$$ which is the factor by which the receiver's clock is faster than the transmitter's. Of course, this is only one of many possible approaches.


Does any error correction code still work in such situation? I'm looking for a kind of error correction code or solution that can correct my codeword in this case:

My message holds $k$ bits, and $2k$ bits codeword (rate is $1/2$) is produced by the generator matrix, e.g. $k=5$, so it yields $10$ bits codeword $\rm 1001100111$. If the codeword is erased for some reason and it becomes $\rm 1xx1x00xx1$ (similar to the codeword transmitted through Binary Erasure Channel, and $\rm x$ is unknown, either '$0$' or '$1$'), can I still correct all the "$\rm x$ bits" in such situation as bit error rate is $1/2$?

I have read about some error correction code, such as block code, convolutional code and LDPC and got the following questions:

All these ECC have pratical application in communication, so do they still work when BER is close to $1/2$ (I think BER is impossible so high)?
Are there any feasible solution to my case? In fact, my message is $30$ bits and I can introduce another $30$ bits redundancy as parity bits. Can I still correct my codeword even if its half bits are erased where the recipient knows which bits he didn't receive properly?
What's the limit of ECC? If it can't correct half bits of codeword that are erased, how many bits can it correct at most?

Any guides or suggestions are appreciated!

I'm afraid that this will be an incomplete answer, but I'm giving it in the hopes that it is better than nothing.

It is impossible to recover any information with a BER of 50% because, as information theory tells us, no information is getting through such a channel.  You can see this intuitively if you realize that no matter what bit you send, 0 or 1, it has equal chance of arriving at the receiver as a 0 or 1, so what you receive tells you nothing about what was actually sent.  In such a case it is theoretically impossible to create an ECC that will help at all.

Erasures are, in some situations, better than errors though.  Reed-Solomon decoders often use erasures to increase their error correction capability.  The reason why it works is because with RS codes you have two sets of unknowns: the error locations and the error magnitudes (the number of errors is also an implicit unknown).  Erasures give the decoder the error location so it only needs to figure out the magnitude.

Talking about "erasures" in the context of bits, though, makes little sense because with bit-based codes there is only one set of unknowns: the error locations.  The error magnitudes are always 1.  So if you have a bit "erasure" you are trading one unknown (error location) for another (error magnitude).

With bit-based codes like Viterbi codes the analog to erasures is using soft-decisions.  Soft-decisions contain "confidence" information- i.e. how sure the receiver is that the bit is really a 1 or really a 0.  That information can be used to greatly improve decoder performance.

All that being said, though, if your error/erasure rate is really that high you will probably not get good performance out of any ECC.


Signal and system convolution In the convolution equation $y(n)=x(n)\circledast  h(n)$ if two things are given like $y(n)$ and $x(n)$ then how can find the $h(n)$?

This is known as deconvolution. A typical approach would be to apply a Fourier transform to both sides
$$
y = x \circledast h \overset{\mathscr{F}}{\Leftrightarrow} Y=XH
$$
and divide both sides by the known $X$ to obtain $H=Y/X$. Then apply an inverse Fourier transform to both sides, yielding $h$.

The linked wikipedia article goes into more detail on variations of this approach.


Calculate interbeat intervals I'm currently reading a paper and I can't seem to make sense of a certain part. A link to the paper: removed

The part I don't quite get is on page 286 (after the part where they explain their filters). It starts with the sentence "Two common methods can be used ...". They explain how they calculate the interbeat intervals of a signal.

The things I do not understand:

  Cubic spline interpolation is applied to the data for ensuring the
  existence of continuous ﬁrst-order differential and second-order
  differential

1. Do they mean that the first & second order derivative of two sequential splines should be the same to ensure that they are smooth ?

Then:

  $m$ is the number of feature points

2. Is this the $n$ in the formula ? ( really none of the symbols are explained of the formula. I doubt t is the value at a certain frame, since the values are really small, which would make it impossible to be in the defined interval.)

Subsequently then:

  We consider the interval between [0.25, 2]

3. what is this interval ? With respect to which variable ?

And this sentence I can't figure out entirely.

  We choose the ﬁrst-order differential point in each cycle of BVP as
  characteristic point to demarcate the interval. The number of
  ﬁrst-order differential points in a period of wave ﬂuctuates narrowly
  than the number of usual maximum points in the presence of measuring
  errors.

4. How do they know anything about any cycle at this point, this is what they are trying to determine at this point ...

5. What is a first-order differential point in a cycle ?

6. Is a feature point a first-order differential point ?

I am trying to implement this (in python), based on this explanation I can't seem to figure out the details to do this.

I splitted the question into smaller bits to allow more precise answers.

Question 1: you seem to not understand what spline interpolation is. The quoted sentence means that you need to use an interpolation function that guarantees the continuity of 1st and 2nd order derivatives. This is the case of cubic spline and a justification for their use.
Question 2: yes, I guess that there is a typo here and $n$ is $m$.
Question 3: it is an interval of sampling rates. It is not very clear though...
Questions 4/5: the paper is clearly not clear or informative at this point. What the authors do is to detect a heart cycle by looking at the 1st order differential points (not defined, but that are opposed to the number of local extrema). My guess is that these points are the zero-crossing of the signal derivative (computed from the spline interpolation, cf. question 1).

As an additional note, the paper is not really well written or informative. For example, Laplacian eigenmaps are poorly exposed, and prior knowledge is required to understand it. Maybe you should pick another reference paper in order to get started in this field.


IQ Sampling and Oversampling Differences? *Confused* Two different thought process that produce two different alias band calculations .... What concepts are these?

If I have two 1 GHz ADC's and they sample with 90 degree offset (A comment said this is Over Sampling)

Each ADC has an Fs/2 of 0.5 GHz
The combined system sample streams have an Fs/2 of 1 GHz because the second ADC captures additional sample points between the samples of the first - doubling the highest frequency possible by one alone.
The system can be thought of as having Fs of 2 GHz, whereas a single equipment Fs is 1 GHz.
After an FFT, the entire spectrum repeats every 2 GHz, the 'system' Fs not the equipment Fs of 1 GHz.
After an FFT, The positive frequencies are 0 to 1 GHz and the negative -1 to 0 GHz.

I am still a bit confused by the thought that maybe everything above is false because maybe IQ sampling means ...

The FFT is being fed a complex number via IQ sample streams.
The FFT fed a complex number does not have negative frequencies - it can span the entire single equipment Fs of 1 GHz.
So Fs/2 is still 0.5 and Fs is still 1, each one is not doubled like the first bullet points suggest because the complex notion of FFT allows frequencies up to the single equipment Fs of 1 GHz when complex sampled.
So also, Fs is 1 GHz and not 2 GHz (3rd bullet point from first set of points) and so the spectrum repeats every 1 GHz and not 2 GHz.
Also in this case, there is no 'spectral inversion' from negative frequencies.

It makes a big difference to where I calculated the alias band and thus my analog filter whether the first set of statements is true or second set is true

Please can you offer advice. Still require assistance if you can 02/03/2020

I think you're confusing IQ sampling with oversampling.  Your first scenario is oversampling, or taking twice the number of samples on a real signal.  Other than stating the sampling takes place 90 degrees apart (it would actually be 180 degrees), your analysis is sound.

In the second section, you describe IQ sampling.  This is usually done by mixing a signal with a nonzero center frequency to baseband using two mixers, which multiply the same signal by a sine and a cosine, respectively.  The two mixing signals, separated by 90 degrees, are generally output from the same oscillator.  The outputs of the two mixers are separately sampled to get the I and Q components of the sampled stream.  This is typically represented as a stream at the original sample frequency but with two values (I and Q) for each sample.

IQ sampling essentially allows one to differentiate between positive and negative frequencies in the baseband signal, so it does in fact double the useful bandwidth of the signal just as oversampling does, but in a different way.  Instead of doubling the Nyquist rate, it revitalizes the band from -Fn to 0.


Does one ever store audio in mp3s with one channel being inverted? As far as I am aware this is something one does for analog processing, not in digital, and one only does it during the actual processing and not in storage (for example in an mp3).

Is my thinking correct?

Well, Audacity has an option to do that on digital files, so I assume that means it is done.

From that page:

  Invert flips the audio samples upside-down, reversing their polarity. The positive samples (above the horizontal zero line in the Audacity Waveform) are moved below the zero line (so becoming negative), and negative samples are made positive.
  Invert does not usually affect the sound of the audio at all, but it can be used for audio cancellation. If Invert is applied to one track and that track is mixed with another uninverted track that has identical audio, the identical audio is cancelled out (silenced).


How to overlap two RAW audio signal I hope I can find some help here.

I have started working on a little project concerning voice recordings. What I would like to is: given two audio files in RAW format (or WAV, as long as I can convert them to arrays of bytes), I would like to overlap them, that it get the RAW signal that is roughly equivalent to playing them at the same time.

Being the complete newbie that I am at signal processing, the first thing I've tried, as soon as I was able to turn my files into arrays, was to average them arithmetically ((x+y)/2) and geometrically (sqrt(x*y)). Clearly both approaces return a shriek from hell. Now I was looking into convolutions, but I don't think it is the way to go. Last I thought I could go with FFT, sum themin the transform space and then convert back, but I'm not quite sure. So here goes:

How can I overlap two RAW audio signals? That is, given two audio signals that sound good by themselves, how can I produce a new audio signal that sounds as if the former were being played at the same time?

Hope it is possible to do it using RAW audio, but other suggestions are appreciated. Thanks

Adding the signals in time domain(sample by sample addition) while ensuring they don't add up to saturate will give you the desired effect. operation is adding and normalizing(like arithmetic mean). in question it is mentioned that taking mean is not giving desired results, this puzzles me. if i understand the desired effect properly then, adding element by element in time domain and normalizing it should do.


PAPR in the SC-FDMA system in wireless communication,compared to the OFDM system , why does the SC-FDMA (single carrier FDMA) system has a lower PAPR (peak to average power ratio)?

regards,
phani tej

This can be explained with the help of some knowledge of probability theory.
Since the output signal is a hybrid signal consists of many(say,Nfft=2048) single carrier signals, which are QAM-modulated by the complex symbols derived from the information bits. Thus leads to(according to the Central Limit Theorem in probability theory) a near Gaussian distributed signal.

As is known to us, Gaussian distribution is not bounded in its range, hence resulting in a high PAPR in OFDM signal. While, in SC_FDMA, things are different. The output signal is essentially a single carrier modulated by the complex symbol, no summation of Nfft signals here. Thus a low PAPR is obtained.

This can also explained intuitively in the time domain. Say, you have a single carrier(in SC-FDMA) modulated by A QAM modulated scheme. The peak of the signal waveform is almost bounded by the carrier amplitude it self. Whilst, in OFDM, where you have a large number of carriers, each with a different frequency. The peak of the output signal is not necessarily bounded by the highest amplitude, since at some time point, all the peaks of the carriers maybe coincide together and thus sum up to form a ultra high peak. While, at some time point, its negative counterpart can be appeared. Although both of the positive ultra high peaks and negative of it are rare in time, it indeed leads to a high PAPR.


Inverse discrete time Fourier transform with differentiation Consider a signal x[n] and its DTFT $X(e^{jω})$ . Assume $X(e^{jω})$ is differentiable. Compute the inverse DTFT of
$j\frac{dX(e^{jω})}{d\omega}$
You should write your answer in terms of $x[n]$ and elementary functions and constants, for example $π/2x[n]$ would be written :
$\pi / 2x[n]$

Hint:
differentiate from Fourier Transform equation:
$$X(e^{j\omega})= \sum\limits_{n=- \infty }^{\infty}x[n]e^{-jn\omega} \Rightarrow \frac{d}{d\omega}X(e^{j\omega})= \frac{d}{d\omega}\sum\limits_{n=- \infty }^{\infty}x[n]e^{-jn\omega}$$
and then find what you want.


Invertibility of Time-Dependent Fourier Transform I am reading Oppenheim & Schafer's (O&S) Discrete Time Signal Processing (2nd or 3rd edition, does not matter) and I find hard to understand a technicality behind the Time-Dependent Fourier Transform (TDFT from now on). Specifically, I have a question regarding the invertibility of the TDFT.

According to the textbook, the TDFT is given by
$$X[n, \lambda) = \sum_{m=-\infty}^{+\infty}x[n+m]w[m]e^{-j\lambda m}$$
which can be thought as the DTFT of a windowed signal portion around sample $n$. So, by moving one sample at a time, we have "a collection" of DTFTs of windowed signal portions around sample $n$, and the TDFT is a 3D-function (consider magnitude only) with axis $n$ vs axis $\lambda$.

Then, O&S mentions that the TDFT is invertible iff $w[m]$ has at least one nonzero sample.

If we take the Fourier synthesis equation, we get: $$x[n+m]w[m]=\frac{1}{2\pi}\int_0^{2\pi} X[n,\lambda)e^{j\lambda m}d\lambda, \: \: \: -\infty < m < +\infty$$
which makes perfect sense as it reconstructs the windowed signal around sample $n$ from complex exponentials weighted by the DTFT $X[n, \lambda)$.
Dividing by the window: $$x[n+m] = \frac{1}{2\pi w[m]}\int_0^{2\pi} X[n, \lambda)d\lambda$$ if $w[m] \neq 0$.

Here's my question: where did the $e^{j\lambda m}$ term go??

Isn't this equation supposed to recover the full windowed signal $x[n+m]$? The integral looks as a number to me (summing $X[n,\lambda)$ values for all $\lambda$ at index $n$ of the TDFT).

The same equation is in all O&S versions, so I guess it's not a typo and I am missing something very fundamental here :D

Thanks in advance.

I have the 2nd edition, and there the case $m=0$ is considered, so the equation reads (Eq. $(10.21)$ on p. 716)

$$x[n] = \frac{1}{2\pi w[0]}\int_0^{2\pi} X[n, \lambda)d\lambda$$

(which of course assumes $w[0]\neq 0$). This makes sense because $e^{j\lambda m}$ equals $1$ for $m=0$.


Using the windowing Method I'm trying to write a C program that can get the PSD of an recorded audio file, I'm using for the fftw3 library. I wrote a program that calculate the PSD so far only for sinus function, but the result is correct, my problem is that I have some trouble to use the windowing method on my data.

Let's say that the audio file has following properties :

fs = 44100;
length = 10sec;

So I got 441000samples, and I want to use windows, doesn't matter which one, with the size of 512, 25% overlapping, and shift my windows for the data. How could I know which value of the result is representing which frequency?

In Matlab it was easy here what I've done for example  :

  [x,fs] = audioread('fileName.wav');
  [xPSD,f] = pwelch(x,hanning(8192),4096,8192*4 ,fs);

Thanks for any help or explanation!

The $k$-th FFT bin will contain the complex amplitude for frequency $\frac{k}{512} \times 44100$.


system response: time vs frequency. Why do I get different magnitudes? I am simulating a linear system described by the following equation

$m \ddot{x} + c \dot{x} + k x = A \sin(\omega t)$

Fs = 1000;      % sampling frequency
tspan = 0:1/Fs:10;
m = 25;
c = 15;
k = 330000;
A = 100;           % forcing amplitude
omega = 114.89;    % forcing frequency

% State space function
odefun = @(t,u)[u(2); 1/m*(-k*u(1)-c*u(2)+A*sin(omega*t))];

[time,u] = ode45(odefun,tspan,[0 0]);

spost = u(:,1);    % displacement
vel = u(:,2);      % velocity

Then I compute the FFT of the displacement:

L = length(spost);
NFFT = 2^nextpow2(L);
SPOST = fft(spost,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);

And this is what I get:

As you can see the maximum value of the displacement is different in the 2 cases, but this does not look right to me. In fact, if I run a structural analysis in the time or in the frequency domain I would expect to obtain the same maximum displacement (since the excitation frequency is equal to the resonance frequency).

How can this be explained?

EDIT

I forgot to say I tried also increasing the simulation time, but still the 2 values do not match

What you are trying to do is to compare the amplitude of your signal with an amplitude of FFT peak, expecting it is the same. The thing is that amplitude of your signal is varying and FFT is taking the "average" over whole duration.

In such case you should use time-frequency representation of your signal, such us spectrogram. It will allow you to observe how amplitude of particular frequency is changing over time. Since you work in MATLAB you can try following code:

win = hamming(256)
noverlap = 0
nfft = 512
X = abs(spectrogram(spost, win, noverlap, nfft))
% Scale by window and multiply by 2 to restore energy from negative freqs
X = 2*X/sum(win)
surf(X)

It will produce something like:

The maximum value of displacement is around $0.053$. Results can be tweaked even further by tweaking even more and using other methods.


SD of a function of signals Say $a_x$, $a_y$, $b_x$, $b_y$ are signals with values in $\mathbb{R}$. Say they are independednt, normally distributed variables with given SDs $\sigma_{a_x}$, $\sigma_{a_y}$, $\sigma_{b_x}$, $\sigma_{b_y}$.

Say $f: \mathbb{R}^4\longrightarrow\mathbb{R}$ is a continuous function of these signals (i.e. $f(a_x, a_y, b_x, b_y)\in \mathbb{R}$ for all $(a_x, a_y, b_x, b_y)\in\mathbb{R}^4$.

General question

What is the SD of $f(a_x, a_y, b_x, b_y)$?
Example/Question

Say $a$ and $b$ are vectors (say in $\mathbb{R}^2$) and $\sphericalangle$ is the angle between $a$ and $b$, say defined as
$$\sphericalangle(a,b) := \mathrm{arcos}\left(\frac{ab}{\vert a\rvert \lvert b\rvert}\right)\quad\text{where}\quad a=(a_x, a_y)\quad\text{and}\quad b=(b_x, b_y).$$

Say $a_x, a_y, b_x, b_y$ are independent and normally distributed, have SDs $\sigma_{a_x}$, $\sigma_{a_y}$, $\sigma_{b_x}$, $\sigma_{b_y}$ respectively.

Then Gauss error proparation yields
$$\sigma_\sphericalangle:= \sqrt{ { \left(\frac{\delta \sphericalangle}{\delta a_x} (a,b) |\mu \ \cdot \sigma_{a_x}\right) }^2 + { \left(\frac{\delta \sphericalangle}{\delta a_y} (a,b) |\mu \ \cdot \sigma_{a_y}\right) }^2 + { \left(\frac{\delta \sphericalangle}{\delta b_x} (a,b) |\mu \ \cdot \sigma_{b_x}\right) }^2 + { \left(\frac{\delta \sphericalangle}{\delta b_y} (a,b) |\mu \ \cdot \sigma_{b_y}\right) }^2}$$

as the error of $\sphericalangle$ (where $..|\mu$ shall note evaluation at average $\mu$).

Question

Is that $\sigma_\sphericalangle$ the SD of $\sphericalangle$?


GMSK Big error in demodulation MatLAB Simulink I am trying to do a PLC simulation system, but my problem is in GMSK demodulator, in my simulation in Simulink MatLAB my output signal is very different from the input signal, below the model:

Below the block parameters:

I start doing this simulation yesterday, the values of each block is the default, and even without de AWGN block the error still appears. Someone knows what is happening?

That’s is because of the delay in the GMSK demodulator. MATLab GMSK demodulator uses Virtebi algorithm. You should delay the Rx BER in order to get back the correct data alignment.


How is the ROC of a transform function determined? Suppose I have $x[n]$ and $y[n]$, and I calculate their respective Z-transforms $X(z)$ and $Y(z)$ as well as their respective ROCs.

Calculating $H(z)$ is as simple as calculating the quotient of $\frac{Y(z)}{X(z)}$, but how do I get $H(z)$'s ROC?

An example:

The ROC of $Y(z)=X(z)H(z)$ is the intersection of the ROCs of $X(z)$ and $H(z)$, unless there are pole-zero cancellations, which change the ROC of the product. In the given example we have to take into account the pole-zero cancellation at $z_0=2$. We can define a new function $X'(z)=X(z)(1-2z^{-1})$ with the pole at $z_0=2$ removed. We do not need to define a new $H'(z)$ with the corresponding zero removed, because that zero doesn't influence the ROC anymore (since the corresponding pole has already been removed). The ROC of $Y(z)$ is now the intersection of the ROCs of $X'(z)$ ($|z|>\frac13$) and $H(z)$. Consequently, the only possible ROC for $H(z)$ is $|z|>\frac23$. This means that the corresponding system is causal and stable, because the ROC is outside the pole with the largest magnitude and the ROC contains the unit circle.


Don't understand coefficients in MATLAB "butter" low pass result I don't unstand the nature of the output in Matlab's "butter" command. Calling

[b,a] = butter(4, 0.5, 'low');

gives me two vectors, which the documentation says are b for numerator, and a for denominator, of the transfer function.

However plotting the frequency response of this transfer function with freqz() does not give me what I expected.

If I plot

freqz(b./a)

I see:

This is not a low pass filter, this is some sort of band stop filter suppressing the middle frequencies. However if I plot

freqz(b)

I see what looks like the low pass filter I was hoping for:

So if I want to use a Butterworth filter to low pass my data, do I only use b? And if so what is the point of a?

You are using wrong syntax. You should call freqz function as:

freqz(b,a)

You do not divide any coefficients. In case of non-recursive filters (FIR) you use only coefficients stored in b vector. When you are using the recursive LCCDE (Linear Constant Coefficients Difference Equation) (mostly IIR filters) then you need also the a vector. I suggest you to set a to 1 if you use non-recursive FIR filters.

You should also watch out for numerical problems with transfer function representation as mentioned in final paragraph of this manual entry and this answer. It's always safer to use z-p-k representation.


Iterative Blind Sinus Signal Suppression There are two real signals in the form of $A_i sin(wt+p_i), i=1,2$. Suppose frequency $w$ of both the signals is the same and amplitude $A_i$ and phase $p_i$ are different. The first signal has unknown $A$ and $p$, but they are constant in the time of observation (for simplicity only, actually they aren't). The second signal is under control. We only have it's difference $x(t) = s_1(t)-s_2(t)$. Some iterative algorithm is needed to minimize power of $x(t)$. So $x(t)$ is available input and $A_2(t)$, $p_2(t)$ are parameters under control.

The question is: is it some robust and rapid algorithm for jointly control both amplitude and phase to minimize power of $x(t)$? My idea was to implement complex gradient method: phase and amplitude at each step are represented as a complex number $C(t)$. Signal $s_2(t)$ can be set as $Re[C(t) exp(iwt)]$. But I can't understand how to define cost function and gradient in this case and also this approach seems to be computationally redundant. Any ideas?

Thanks in advance!

I assume you do everything in discrete time. If I understood correctly then the observed signal is

$$x(n)=s(n)-A\sin(\omega_0n+\phi)$$

where $s(n)$ is a sinusoidal signal with frequency $\omega_0$ and unknown amplitude and phase. The parameters $A$ and $\phi$ are to be optimized such that the power of $x(n)$ is minimized. In a real-time algorithm you can define the cost function as the sum of $x^2(n)$ over a fixed length $N$:

$$F(A,\phi,n)=\sum_{k=n-N+1}^{n}x^2(k)$$

where $n$ is the current sample. For a steepest descent algorithm you need the gradient:

$$\frac{\partial F}{\partial A}=-2\sum_{k=n-N+1}^{n}x(k)\sin(\omega_0k+\phi)\\
\frac{\partial F}{\partial \phi}=-2A\sum_{k=n-N+1}^{n}x(k)\cos(\omega_0k+\phi)$$

The parameters are now updated in the following way:

$$A_n=A_{n-1}-\mu\frac{\partial F_{n-1}}{\partial A}\\
\phi_n=\phi_{n-1}-\mu\frac{\partial F_{n-1}}{\partial \phi}$$

where $\mu$ is the step size. Note that in every iteration you need to recompute the gradients with the current values of $A$ and $\phi$.

The figure below is a simulation of this gradient algorithm with the unknown signal $s(n)=2\sin(\pi n/3+\pi/2)$, $N=50$, $\mu=0.01$, and with zero initialization of both amplitude and phase (amplitude $A$ in blue, phase $\phi$ in green).


Tracking object that is missing in some frames I am trying to track small object in videos using LK
but the problem is, some time the object does not appear in the frames !

How can I track an object that does not appear in few frames ??
 how can I approximate it is place automatically ,,,

I am trying to find the object using Optical flow LK ... but this technique lost when my   object does not appear

Any suggestion ideas ???

The short answer is: use your optical flow algorithm to get the best-fit warping of some number of previous images into your current image.  Then assume that your object is moving on a line through the images, and do some kind of robust best-fit to that line.  When the object reappears it is more likely to be near the predicted projection of the line than not.

You didn't give enough specifics in your question to really say more than that.

I assume by "LK" you mean "the Lucas-Kanade optical flow algorithm?"  But what are you really doing?  Are you taking image pyramids of each image and then doing multiple iterations of warping at each level of the pyramid?  Are you really doing the simplest form of Lucas-Kanade or are you actually doing one of the varieties of Kanade-Lucas-Tomasi where you are concentrating on "corners"?

What is your description of an "object"?  Is it really just a "feature" (a corner point) or a small set of features?

You may need to use a more accurate and robust optical flow algorithm.  For example, some variation that does dense estimation (rather than the sparse estimation done by Lucas-Kanade) and that uses more robust distance metrics than least squares.

For example:

Secrets of optical flow estimation and their principles
Sun, D., Roth, S., and Black, M. J.,
IEEE Conf. on Computer Vision and Pattern Recog., CVPR, June 2010.


Complex numbers in $\tt ifft$ of an MMSE amplitude estimator I am trying to reconstruct a signal from a noisy speech using an MMSE algorithm proposed long time ago by Ephraim and Malah (1984). After going through the algorithm, I got a matrix A which represents the magnitude of the reconstructed signal. With the help of this group, I now know that I need to use this magnitude with the phase of the noisy signal in order to recover the signal. Now the problem is that after doing an ifft, I get complex numbers.

Does this mean there is a mistake in implementing the algorithm?
I read somewhere that I have to have symmetry in my input to an ifft function, how do I ensure such symmetry is applied here?

Here is my MATLAB command and the equation

A = A_hat.*exp(1i*nu);
x_new = ifft(A);

  \begin{align}
\hat A_k&=\Gamma\left(1.5\right)\frac{\sqrt{v_k}}{\gamma_k}M\left(-0.5;1;-v_k\right)R_k\\
&=\Gamma\left(1.5\right)\frac{\sqrt{v_k}}{\gamma_k}\exp\left(-\frac{v_k}{2}\right)\left[\left(1+v_k\right)I_0\left(\frac{v_k}{2}\right)+ v_k I_1\left(\frac{v_k}{2}\right)\right]R_k.\tag{7}
\end{align}

Note that the MMSE estimator gives you a real-valued gain function with which the noisy signal amplitude is modified. Equation (7) is actually

$$\hat{A}_k=G_kR_k\tag{1}$$

where the gain $G_k$ is a function of $v_k$ and $\gamma_k$. You want to reconstruct the original signal as

$$\hat{X}_k=\hat{A}_ke^{j\phi_k}\tag{2}$$

where $\phi_k$ is the phase of the noisy signal $Y_k=R_ke^{j\phi_k}$. Combining (1) and (2) gives

$$\hat{X}_k=G_kR_ke^{j\phi_k}=G_kY_k\tag{3}$$

So you can simply multiply the complex spectral component of the noisy signal with the gain $G_k$ to get the estimation of the spectral component of the input signal.

In Matlab this is easily done. You just have to extend the gain vector to satisfy the symmetry condition for real-valued signals:

X = Y.*[gain; gain(N-1:-1:2];

where I assume that the vector gain has elements $G_k$, and that N is its length. The relation between the FFT length and N is N = Nfft/2+1.


SQNR when multiplying 2 16 bit signed integers and truncating result? I'm trying to figure out what the Signal to Quantized Noise Ratio (SQNR) is if I multiply 2 16-bit signed integer values and store the result as a 16-bit signed integer, dropping the least significant bits. Since this is signed multiplication, I'm dropping 15 bits of precision.

For example:

a (16 bit signed) * b (16 bit signed) = c (16 bit signed) without the LSBs.

And to be more clear about what's after this, I am summing many successive values of the resulting 'c' - this sum is stored in an infinite length signed integer representation.

How do I calculate the SQNR in this situation?

When you multiply the two 16 bit numbers you get the "true" 32 bit number and then you truncate to get a 16 bit number.  The noise that is introduced when you truncate is a random value from $0$ to $2^{16}-1$ ($65535$), which averages out to $32767.5$ and has an average power (average of the square) of $1.43E9$.  If you were rounding instead of truncating it would be a random value from $-(2^{15})$ to $2^{15}-1$, which averages out to approximately $0$ (no bias), and has an average power of $3.58E8$.  That's why it's better to round than truncate, but that's a story for another day.

The truncation noise, therefore, is $1.43E9$.  You get the signal power by multiplying your truncated data by $2^{16}$ (to make the 16 bit values into 32 bit values), squaring them, and then calculating the mean.  You then calculate the SQNR by dividing the signal power by the truncation noise, $1.43E9$.


Basic sampling excercise - unsure of rationale informing np.arange parameters for time array I'm extremely new to DSP, so I apologise if it's a very basic question.

I'm having issues understanding the following exercise:

"Create an array of 201 samples containing a sin wave with a frequency , f, of 2 Hz, i.e. the period is 0.5 seconds. Use a sampling rate, fs of 100 Hz, i.e 100 samples per second. Use the numpy method arange to generate the array. The first step is to declare the variables for f and fs and then create an array of times for the samples in a 2 second range, i.e.

t = np.arange(−1,1+1/fs,1/fs)

which gives a array starting at -1 and ending at 1.01."

I get how you create a range in Python, and I get that the parameters are the starting value, the end value, and the steps. However, I'm not quite sure what this is doing: why does it start at -1? Why does it end at 1+1/fs? Why are the time steps 1/fs? And lastly, what exactly is this giving me? Just the times at which the wave is sampled?

I've tried Googling, but I don't really know what I should be looking for. Even a pointer to the right keyword to Google to answer those questions would be greatly appreciated.

The Hz values are used to get your parameters for defining your signal.

I don't think there is any reason that they went from -1 to 1 instead of, say, 0 to 2.  Notice that the end point is not included in the range.

The spacing in your sampling is 1/100 of a second.  Yes, the array is giving you values on a time scale.

To generated the signal your are going to need:

f = 2
signal = np.sin( f * 2*np.pi * t )

The $2\pi$ comes as a unit conversion of radians per cycle.  Therefore the argument, in units is:

$$ \frac{cycles}{second} \cdot \frac{radians}{cycle} \cdot seconds = radians $$

Hope this helps.

Ced

================================

Followup:

Here is a sample program.

import numpy as np
import matplotlib.pyplot as plt

f = 2.0
fs = 100.0

t = np.arange( -1.0, 1.0+1.0/fs, 1.0/fs )
signal = np.sin( f * 2*np.pi * t )

plt.plot(t, signal, c='g')
plt.show()

Notice that I added a bunch of ".0"s to force Python into using floats.  Without it, the program passed 1/fs as the integer 0 to arange.


How CDMA receiver extract it's corresponding data from the receiving modulated & superposition-ed signal? I have difficulty understanding CDMA.

I have read all those PN/walsh code things, and understand orthogonality very well.

But none of the resources I found demonstrate how a receiver end extract data from modulated & superposition-ed signal.

Most of them just stop at the digital encoding & decoding part and skip the modulation part.

But in real implementation signals of different walsh code are modulated then superposition-ed.

I cannot figure out how it is done.

======================================================

P.S.
I am confused that does the sender add the signals before or after modulation?

If before, take BPSK as an example, there is a chance there would be >2 different values in the added result vector, which can no longer be encoded in BPSK.

If after, then the signals are simply super-positioned, the receiver would often see both 0&1, with different amplitudes, at the same time, i.e. Acos(wt)+Bcos(wt+pi), which seems a bit hard to decode.

You clarified that your question really is how multiple users can be separated in code space. Since you said that you understand orthogonality well, I assume you just need an example to tie it all together.

With reference to the below examples of the DFT (OFDM) and Hadamard Matrices (CDMA) where each row of either represents a different orthogonal code that we can assign to a different user, channel or resource.  For CDMA and with reference to the 8th order Hadamard Matrix, which I am referring to as $H_8$, we have a code set of 8 possible codes which we can use for multiple access with 8 different users by assigning each user a different code.

We can then choose a modulation to use for each symbol, for example BPSK, QPSK, QAM etc. BPSK is simplest so I will start with demonstrating that assuming one particular user that has the 2nd row as their code which is $H_{8,2} =\begin{bmatrix}1 &-1 &1 &-1 &1 &-1 &1 &-1\end{bmatrix}$:

For BPSK the users data such as 0 1 1 0 1 1 1 0 is mapped to symbols -1 and 1 (bi-phase 0 or 180°). To mix with the code we simply multiply each symbol with the code sequence for that user, so to send the data above we would do:

User Data: $\begin{bmatrix}0 &1 &1 &0 &1 &1 &1 & 0\end{bmatrix}$

Data Symbols: $\begin{bmatrix}-1 &1 &1 &-1 &1 &1 &1 &-1\end{bmatrix}$

Symbol 1: $-H_{8,2} =\begin{bmatrix}-1 &1 &-1 &1 &- 1 &1 & -1 &1\end{bmatrix}$

Symbol 2: $H_{8,2} =\begin{bmatrix}1 &-1 &1 &-1 &1 &-1 &1 &-1\end{bmatrix}$

Symbol 3: $H_{8,2} =\begin{bmatrix}1 &-1 &1 &-1 &1 &-1 &1 &-1\end{bmatrix}$

Symbol 4: $-H_{8,2} =\begin{bmatrix}-1 &1 &-1 &1 &- 1 &1 & -1 &1\end{bmatrix}$

Symbol 5: $H_{8,2} =\begin{bmatrix}1 &-1 &1 &-1 &1 &-1 &1 &-1\end{bmatrix}$

Symbol 6: $H_{8,2} =\begin{bmatrix}1 &-1 &1 &-1 &1 &-1 &1 &-1\end{bmatrix}$

Symbol 7: $H_{8,2} =\begin{bmatrix}1 &-1 &1 &-1 &1 &-1 &1 &-1\end{bmatrix}$

Symbol 8: $-H_{8,2} =\begin{bmatrix}-1 &1 &-1 &1 &- 1 &1 & -1 &1\end{bmatrix}$

Thus the users 7 data bits in this example was converted to 56 actual bi-phase symbols to be transmitted (in this case "spread spectrum" specifically since the data rate and hence the occupied bandwidth was increased). To receive the sequence for this user, we simply correlate (multiply and accumulate) the received sequence to each code and detect which one is strongest and what the magnitude and phase is to demodulate the underlying modulation.

Since the codes are orthogonal, you could encode a 2nd users data in all the same time slots, add the resulting encoded symbols and we would be able to separate the users in the receiver by correlating to each code in turn.

For example if we want to send a data bit 0 to User 1 who uses the 2nd row of $H_2$ as we did above:

$-H_{8,2} =\begin{bmatrix}-1 &1 &-1 &1 &- 1 &1 & -1 &1\end{bmatrix}$

And at the same time we send a data bit 1 User 2 who uses the 1st row of $H_2$:

$H_{8,1} =\begin{bmatrix}1 &1 &1 &1 &1 &1 & 1 &1\end{bmatrix}$

The combined result that we would transmit would be the addition of the two:

$t = \begin{bmatrix}0 &2 &0 &2 &0 &2 & 0 &2\end{bmatrix}$

In the receiver we can correlate to each user by taking the dot product (multiply and accumulate which is correlation) of the time aligned received sequence with each code:

$t \cdot H_{8_1} = +8$  The best estimate is  +1 symbol was transmitted = data bit 1

$t \cdot H_{8_2} = -8$  The best estimate is a -1 symbol was transmitted = data bit -1

When we don't have time aligned in the receiver, we can use special acquisition sequences that have good autocorrelation properties for finding time offset. This also holds true when we can't guarantee time alignment between all the codes such as when the users are transmitting- Walsh codes are orthogonal when time aligned to each other, but do not have good cross correlation properties when there are time offsets between the codes.

If the modulation of choice was QPSK, then the users data would be mapped 2 bits at a time to 4 symbols, for example $1 0 \rightarrow 1$, $1 1 \rightarrow j$, $0 1 \rightarrow -1$ and $0 0 \rightarrow -j$ but the process would be exactly the same.

User Data: $\begin{bmatrix}0 &1 &1 &0 &1 &1 &1 & 0\end{bmatrix}$

Data Symbols: $\begin{bmatrix}-1 &1 &j &1\end{bmatrix}$

Symbol 1: $-H_{8,2} =\begin{bmatrix}-1 &1 &-1 &1 &- 1 &1 & -1 &1\end{bmatrix}$

Symbol 2: $H_{8,2} =\begin{bmatrix}1 &-1 &1 &-1 &1 &-1 &1 &-1\end{bmatrix}$

Symbol 3: $jH_{8,2} =\begin{bmatrix}j &-j &j &-j &j &-j &j &-j\end{bmatrix}$

Symbol 4: $H_{8,2} =\begin{bmatrix}1 &-1 &1 &-1 &1 &-1 &1 &-1\end{bmatrix}$

We can then multiply the resulting real or complex numbers from the coding above with a carrier frequency to translate the modulated signal to any frequency carrier of interest.

Everything involved to actually receive this properly is more complex due to adjusting for timing offset, frequency offset, phase ambiguity, channel distortion, phase noise, etc etc but that aspect is similar across all types of receivers. With regards to removing the code modulation we simply multiply again by the code once time aligned.

Observe the similarity between OFDM and CDMA in that the DFT is just another orthogonal code space. In the case of Walsh Codes we use two symbols such as 0 and 1, while in the case of the DFT we use the Nth roots of unity which are complex points equally distributed on the unit circle. Both allow us to transmit to multiple users (or multiple channels to one user) in the same time slot, in contrast to TDMA. The DFT is more intuitive since we are so used to visualizing different frequencies such that we see how they don't overlap, but the same thing occurs with Walsh codes it's just that we are in a different code space. With that in mind, compare carefully the 2, 4 and 8 point DFT matrices to Walsh code sets of the same order (Hadamard Matrices).

Often for Walsh Codes, Hadamard Matrices or order N are designated as $W_N$ but to avoid confusion with the DFT twiddle factors $W_N^{nk}$ I will use $H_N$ here as Hadamard Matrices:

2 Pt DFT Matrix

$\begin{bmatrix}1& 1\\1 &-1 \end{bmatrix}$

2 Pt Hadamard Matrix

$H_2 = \begin{bmatrix}1& 1\\1 &-1 \end{bmatrix}$

4 Pt DFT Matrix

$\begin{bmatrix}1& 1 & 1 & 1\\1 & -j & -1 & j \\1& -1 & 1 & -1\\1 & j & -1 & j \end{bmatrix}$

4 Pt Hadamard Matrix

$H_4 =\begin{bmatrix} H_2 & H_2 \\ H_2 & -H_2\end{bmatrix} = \begin{bmatrix}1& 1 & 1 & 1\\1 &-1 & 1 & -1 \\1& 1 & -1 & -1\\1 &-1 & -1 & 1 \end{bmatrix}$

8 Pt DFT Matrix

$\begin{bmatrix}1& 1 & 1 & 1 & 1& 1 & 1 & 1\\1 & W_8^1 & W_8^2 & W_8^3 & -1 & -W_8^1 & -W_8^2 & -W_8^3\\1 & W_8^2 & -1 & -W_8^2 & 1 & W_8^2 & -1 & -W_8^2\\1 & W_8^3 & -W_8^2 & W_8^1 & -1 & -W_8^3 & W_8^2 & -W_8^1 \\1& -1 & 1 & -1 & 1& -1 & 1 & -1\\1 & -W_8^1 & W_8^2 & -W_8^3 & -1 & W_8^1 & -W_8^2 & W_8^3\\1 & -W_8^2 & -1 & W_8^2 & 1 & -W_8^2 & -1 & W_8^2\\1 & -W_8^3 & -W_8^2 & -W_8^1 & -1 & W_8^3 & W_8^2 & W_8^1 \end{bmatrix}$

8 Pt Hadamard Matrix

$H_8 =\begin{bmatrix} H_4 & H_4 \\ H_4 & -H_4\end{bmatrix} =\begin{bmatrix}
1& 1 & 1 & 1 &1& 1 & 1 & 1\\
1 &-1 & 1 & -1 &1 &-1 & 1 & -1  \\
1& 1 & -1 & -1 &1& 1 & -1 & -1\\
1 &-1 & -1 & 1 & 1 &-1 & -1 &1 \\
1& 1 & 1 & 1 &-1& -1 & -1 & -1\\
1 &-1 & 1 & -1 &-1 &1 & -1 & 1  \\
1& 1 & -1 & -1 &-1& -1 & 1 & 1\\
1 &-1 & -1 & 1 & -1 &1 & 1 &-1\end{bmatrix}$

Where in the 8pt DFT matrix, $W_N^{n} = e^{-j2\pi n/N}$

We tend to give a more physical preeminence to the concept of frequency separation as if that is more representative of the real world than code separation, and think of code separation as just a mathematical construct. They are both equivalently mathematical constructs and both will result in waveforms that will occupy the same spectral bandwidth as a transmitted signal (and there we have decided to use frequency space to share our use of the ether overall). Step back for a moment and ignore all your training and experience with sine waves and observe what is happening above and how similar they are, and how they both have the same level of physical meaning if you are trying to associate it with one. This is similar to the unfortunate naming of "real" and "imaginary" numbers--- one is no more "real" than the other! I have had people tell me that they believe complex numbers don't actually exist like real numbers do- you can observe a real number in the lab with a scope probe--- well you can observe a complex number in the lab with two scope probes! It is all equivalently math that we use to describe our physical world.

My prior answer before the question was clarified which I will retain since it is related:

I don't think your question is specific to CDMA but more in general what does mixing and modulating mean. What I am about to say is applicable to CDMA but also to any receiver that extracts data from a "mixed and modulated" signal.

Modulation is specifically the process of mapping data to a symbol for transmission, for example BPSK, QSPK, QAM, FSK, MSK, etc are all modulations. For one of the simplest, BPSK, we map data bits from 0 to 1 to phases of the carrier 0° and 180°.

Mixing usually refers to the process of translating a signal from one frequency to another using time domain multiplication (devices that do this in the analog domain are called "mixers", but we can also do this in the digital domain but typically refer to them simply as multipliers - they are doing the same thing).

To keep this simple I am not going to get into complex signals but will show the process of modulating and mixing a BPSK signal.

Data = 0 1 0 1 1

BPSK Modulation :  -1 +1 -1 +1 +1

Multiply with a carrier frequency ($\omega_c$):

To send =1:   $(-1)\cos(\omega t)$

To send +1:   $(+1)\cos(\omega t)$

To demodulate multiply with local oscillator with frequency ($\omega_c$) that has been locked the carrier using a carrier recovery loop:

If receiving a =1:

$(-1)\cos(\omega t)\cos(\omega t) = \frac{(-1)}{2}\cos(2\omega t)+\frac{(-1)}{2}\cos(0)$

Filter out the high frequency component to get $-\frac{1}{2}$  This maps to the 0 bit received.

If receiving a +1:

$(+1)\cos(\omega t)\cos(\omega t) = \frac{(+1)}{2}\cos(2\omega t)+\frac{(+1)}{2}\cos(0)$

Filter out the high frequency component to get $\frac{1}{2}$ This maps to the 1 bit received.


Finding maximum using DFT I'm trying to find an efficient way to compute maximum of a signal using its DFT. More formally:

$$\max\left\{ \mathcal F^{-1}\left(X_k\right)\right\}, X_k\text{ is the DFT of the signal and } \mathcal F^{-1} \text{ is the IDFT}$$

The original signal $x(n)$ is real, and it has some noise when sometimes there are wide peaks.

I'm looking for a solution that is quicker than the actual IFFT since the signal is very long (but we have its DFT).

This sounds like a good candidate for the Sparse Fast Fourier Transform (sFFT). Have a look at H. Hassanieh, P. Indyk, D. Katabi, and E. Price, sFFT: Sparse Fast Fourier Transform, 2012. I don't know details of the algorithm. You have the domains swapped, which should be fine as the discrete Fourier transform (DFT) and its inverse (IDFT) are nearly identical.

sFFT has been discussed also here on dsp.stackexchange.com.


Why Coherence is not a valid metric when performing impact excitation? INTRODUCTION
I have understood that Coherence is a function that explains the linear relationship between an excitation signal and a response signal. I know how it is calculated and why it is bounded between 0 and 1.
In a document by Bruel&Kjaer (page 38) I found the following:

The deterministic character of impact excitation limits the
use of the Coherence Function**. The Coherence Function will show a "perfect" value of 1 unless:

There is an antiresonance, where the signal-to-noise ratio
is rather poor. No particular attention needs to be
paid to this. Taking a number of averages should make
the FRF curve smooth (for noise at the output choose
H1).

The person conducting the test impacts the structure in
a scattered way, with respect to point and direction. This
should be minimized so that the Coherence is higher
than 0,95 at the resonances. If the impact point is close
to a node point the Coherence may be extremely low
(≈ 0,1). This is acceptable however, since the modal
strength at this point is weak, and not important for the
analysis.

My worry arose since I was getting a value of "1" for all the frequencies when analyzing the coherence for some impact testing evaluation.
QUESTION:
Why Coherence is not useful when evaluating the relationship between input and output, for experiments where impact excitation is used? Why do I always get a "perfect 1"?

a good reference is:

  Carter GC. Coherence and time delay estimation. Proceedings of the IEEE. 1987 Feb;75(2):236-55

The coherence is actually complex valued and bound by a magnitude of one.  I believe you are referring to the magnitude squared coherence as used in Carter's paper.  The term coherence by itself is used in optics.

The problem with the MSC is the probability density of the estimate of MSC, it takes a lot of independent samples to get a reasonably accurate estimate. When applied to a single sample, the confidence region is [0 1] at any type two error, which essentially implies that a single sample has to considered a deterministic signal and not the random model.  The following paper covers how to calculate the pdf of the sample MSC.

  T. Barnard, "Legendre polynomial expressions for the probability
  density function of magnitude-squared coherence estimates," in IEEE
  Transactions on Acoustics, Speech, and Signal Processing, vol. 29, no.
  1, pp. 107-108, February 1981. doi: 10.1109/TASSP.1981.1163516
  Abstract: Carter, Knapp, and Nuttall [1] presented a discrete Fourier
  transform (DFT) method for estimating the magnitude-squared coherence
  between two zero-mean wide-sense-stationary random processes. This
  paper shows that a simple Legendre polynomial expression can replace
  the hypergeometric function used in the probability density function
  expression of Carter et al. [1] (see also Fisher [2] and Goodman [3])
  for the coherence estimate obtained from nonoverlapped DFT's. A
  standard Legendre polynomial recursion provides a recursive method for
  evaluating the density expression. keywords: {Polynomials;Probability
  density function;Density functional theory;Frequency;Discrete Fourier
  transforms;Signal processing algorithms;Random
  processes;Instruments;Detectors;Distributed computing}, URL:
  http://ieeexplore.ieee.org.mutex.gmu.edu/stamp/stamp.jsp?tp=&arnumber=1163516&isnumber=26151

If you have the symbolic toolbox, you can use the following matlab code to calculate the MSC pdf using the expression from the Carter paper.

clear all
%nd=128
n=[ 1 32 512 768 ]
c=linspace(0,1,512);
p=zeros(length(n),length(c));
for k=1:length(n)
nd=n(k)
%syms  C Ch
C=.5;
p(k,:)=Magsquarecoherepdf(c,C,nd);
end
figure(1)
plot(c,p,'Linewidth',2)
h1=legend('1', '32','512', '768');
title(h1,'Number Trials');
title(['Sample Density given $\mid \gamma \mid^2 = $',num2str(C)],'FontWeight','bold','Interpreter','latex')
h=xlabel('Sample $\mid \tilde{\gamma} \mid^2$','FontWeight','bold','Interpreter','latex');

function [p]=Magsquarecoherepdf(c,C,nd)
%
p=zeros(size(c));
for i=1:length(c)
Ch=c(i);
t1=log10(hypergeom([1-nd 1-nd],1 , Ch*C));
assert(~isnan(t1),'t1 nan')
t2=log10((nd-1).*((1-Ch.*C)./(1-C).^2));
assert(~isnan(t2),'t2 nan');
t3=log10(((1-C).*(1-Ch)./(1-Ch.*C).^2))*nd;
assert(~isnan(t3),'t3 nan');
p(i)=10.^(t1+t2+t3);
%assert(~isnan(p(i)),['p  assert ',num2str(t1),' ', num2str(t2),' ', num2str(t3)])
end
return

If you repeat your impacts, you have an ensemble and I believe you could apply the MSC to your problem.

edit:

If you might permit an analogy. lets say you are asked to determine if a coin is fair in a toss. You can't say very much from a single flip of a coin. you have to repeat the toss , and keep on repeat tossing it. the number  of tosses required for a specific accuracy is determined from a probability law. a single toss tells you nothing.

the same is true for the MSC. a single trial says nothing. it needs to be repeated and the probability density in the code above specifies the number of trials required for  a specified accuracy


Equiripple Filter-Design I have a requirement to design minimax filter with linear programing (linprog in MATLAB).
To build the filter I must choose vector $\omega$ (frequency sample), how many values I need to take to get the optimal result? How to spread them across the interval $[0,\pi]$?

There is no real rule, but you would usually choose around $10N$ frequency points, where $N$ is the desired filter order. Distribute the points equidistantly over pass band(s) and stop band(s). This does not mean that they are equidistant in the interval $[0,\pi]$ because usually there is at least one "don't care" region where no desired response is specified (such as in the transition band between pass band and stop band).


Possibilities to get the amplitude of a signal at specific frequency from IQ I'm newbie to the dsp topic and my understanding of math symbols is limited.
So at the moment I'm little bit lost to find the right resources.

I receive an IQ stream from a SDR.
Sample rate: 250k/s
The signal I'm interested in is at -10kHz.
I need the amplitude of only the -10kHz frequency.

I've found out how to shift the IQ data frequency up/down. But I thought -10kHz is a nice place to be away from DC, and for buffer sizes and sample rate.

What is the most simple solution? Are there other possibilities next to one of the Fourier transformations? Looks like overkill for my simple signal.

But I've found lot of theory stuff handling IQ signals with only one frequency, but nothing enlightening to understand the demodulation of one frequency out of a mix of many from IQ.

I know about the phase change. But something I miss to understand how to collect and calc the right data for my frequency from the IQ.

What would be a nice book (German/English) which explains the stuff without getting lost in math vocab and theory after a few pages? Anything from a software developers perspective?

Thx.

The easiest approach is to compute the one DFT bin of interest directly, which is a correlation to the frequency of interest by complex conjugate product and sum. Given you stated -10KHz frequency and with IQ samples then this frequency can be described using the exponential:

$$f_{in} = e^{-j\omega_n t} = I-jQ$$

Where

$I = \cos(\omega_n t)$

$Q = \sin(\omega_n t)$

$\omega_n = 2\pi10 \text{ KHz}$

Expressing this in samples with $f_s = 250 \text{ KHz}$:

$I_n = \cos(\omega_n/f_s n) = \cos(2\pi n/25) $

$Q_n = \sin(\omega_n/f_s n) = \sin(2\pi n/25) $

The normalized correlation would give you the magnitude and phase:

$$X[\omega_n] = \frac{1}{N}\sum_{n=0}^{N-1}x[n]e^{-j\omega_n} = \frac{1}{N}\sum_{n=0}^{N-1}x[n]\cos(2\pi n/25) + j\frac{1}{N}\sum_{n=0}^{N-1}x[n]\sin(2\pi n/25) $$

or equivalently in terms of the complex samples $I+jQ$:

$$X[\omega_n = -2\pi n/25] = \frac{1}{N}\sum_{n=0}^{N-1}x[n]e^{j2\pi n/25}$$

Using an integer multiple of samples will result in an exact result as it would be the equivalent of a bin center for the DFT:

$$ N = m\frac{250 \text{ KHz} }{10 \text{ KHz} } = 25m$$

for any integer $m>0$

The FFT performs this for all bins with $N\log_2(N)$ operations while here $N$ operations are given, thus when only one bin is needed computing that bin using the DFT directly results in a savings of $log_2(N)$. When multiple but still relatively small number of bins are needed, the Goertzel algorithm is a good choice.


Fourier transform of triangular function Determine $X(\omega)$.

$g(t)$: I understand how to create a box from [-1,1] of amplitude 1/2.
$x(t) = g(t) * g(t)$
$X(\omega) = G(\omega)G(\omega)$

the solution I am seeing says that $G(\omega) = \frac{2\sin(\omega)}  {2\omega}$

I don't understand where $\sin$ came from and that the values of the 2s correlate to. I've seen proofs, but can someone provide a simple explanation of what the variables are. Thanks

A triangular function can be generated by convolving two box functions as shown below.

This is where your Step 2 comes from.

The fourier transform of a convolution $g(t) \ast g(t)$ can be calculated by multiplying the fourier transform of $g(t)$ with itself, i.e. $G(\omega)G(\omega)$.

Recall that the Fourier transform of a box function is a Sinc function ($\textrm{sinc}(x) = \frac{\textrm{sin}(x)}{x}$).

Hence, $G(w)$ is some scaled version of a sinc function, and the Fourier transform of the triangular function is $G(w)^2$.


Integer upsampling - leave original samples unaltered I'm looking at the problem of integer upsampling (i.e. interpolation). This can be achieved by first increasing the sample rate by factor $L$ (i.e. inserting $L-1$ zeros between the original samples) and then low-pass filtering the signal (see e.g. Oppenheim/Schafer, "Discrete Time Signal Processing", Sect. 4.6.2).

When I compare the the samples of the original signal with the corresponding samples of the upsampled signal, I observe a slight difference in amplitude. Is there a way to avoid this and to leave the original samples unaltered? That's what I would expect of an interpolation method.

The filter I used is a symmetric FIR filter with an even number of coefficients and for comparison I compensate the delay introduced by that filter.

I know that there are other interpolation methods (linear, cubic, hermite, etc.), but I would like to use this particular method.

I hope I could make myself understood. Thanks in advance for your help!

You need to be careful with the delay compensation. If your symmetric filter has an even number of taps, then you get a non-integer delay. If you use a symmetric filter with an odd number of coefficients then the delay of $\frac{N-1}{2}$ samples is an integer number ($N$ is the number of filter coefficients).


Flat fading criterion of OFDM subcarrier spacing I have always thought that the OFDM subcarrier spacing $\Delta f$ is chosen such that

not too small because Doppler spread can destroy subcarrier orthogonality
not too large to avoid Cyclic Prefix (CP) overhead because OFDM symbol period $T_u = 1/\Delta f$ needs to be much larger than $T_{CP}$ which, to avoid OFDM symbol ISI, must be larger than multipath delay spread $\tau_m$ which depends only on the given physical environment.

The second criterion implies that $T_u = 1/\Delta f \gg T_{CP} > \tau_m \implies \Delta f \ll 1/\tau_m \sim B_c$ where $B_c$ is coherence bandwidth, this is flat fading definition. I means that flat fading is a consequence of the design criterion, not a criterion itself.

But in this wikipedia article Wikipedia's Fading Article, it is said that

  Since different frequency components of the signal are affected independently, it is highly unlikely that all parts of the signal will be simultaneously affected by a deep fade.Certain modulation schemes such as orthogonal frequency-division
  multiplexing (OFDM) and code division multiple access (CDMA) are
  well-suited to employing frequency diversity to provide robustness to
  fading. OFDM divides the wideband signal into many slowly modulated
  narrowband subcarriers, each exposed to flat fading rather than
  frequency selective fading.

It seems that flat fading is the design objective.

If a portion $B_c$-wide is in deep fade, all subcarriers in that $B_c$ are in deep fade, I am not sure that increasing diversity is the purpose. I mean : could  50 errors per 100 bits be better for channel coding than 5 error per 10 bits?

I understand that $\Delta f < B_c$ facilitates channel estimation in frequency domain by adding at least one pilot per $B_c$. Thus my question : is channel estimation the only reason to choose $\Delta f$ as small as possible ?

Now I look at the time domain, $\Delta f > B_c \implies T_u < \tau_m$. As being visualized in figure below, nothing bad happens. We can capture all the energy from all physical paths. Is there any explaination if we look only at the time domain?

Update : I am doing some math to see how the flat fading condition leads to single-tap equalizer.

Let $N$ the number of subcarrier, or FFT size (no zero subcarrier)
$T_s$ is sampling period $\implies N\Delta f = 1/T_s$ or $\Delta f T_s = 1/N$

The signal at TX after IFFT :

$$x(t) = \frac{1}{N}\sum\limits^{N-1}_{k=0} x[k]e^{j2\pi k \Delta f (t - N_{CP}T_s)}\text{ with }0 \leq t < (N_{CP} + N)$$

The received signal from $L$ multipath propagation:

$$y(nT_s) = \sum\limits^{L-1}_{l=0}h_l(n T_s)x(n T_s - l T_s)\text{ with }0 \leq n \leq N_{CP} + N-1$$

Take the part of FFT window as in the figure above, $N_{CP} \leq n \leq N_{CP} + N-1$. Set $m = n - N_{CP}$ and suppose that $h_l(t)$ is time-invariant in $0 \leq t < (N_{CP} + N)$ so that the time index of $h_l(t)$ can be dropped (this is the underspread assumption which is valid for typical channels) :

$$\begin{align}
y[m] &= y(m T_s) \\
&= \sum\limits^{L-1}_{l=0}h_l \sum\limits^{N-1}_{k=0} x[k]e^{j2\pi k \Delta f (m - l)T_s}\\
y[m] &= \sum\limits^{L-1}_{l=0}h_l \frac{1}{N}\sum\limits^{N-1}_{k=0} x[k]e^{j2\pi k (m - l) / N}
\end{align}$$

with $0 \leq m \leq N - 1$

Note that $T_s\Delta f = 1/N$.

Take FFT of $y[m]$:

$$\begin{align}
z[k_0] &= \sum\limits^{N-1}_{m=0} y[m] e^{-j2\pi k_0 m / N} \\
&= \sum\limits^{N-1}_{k=0}x[k] \sum\limits^{L-1}_{l=0} h_l e^{-j2\pi k l/N} \frac{1}{N} \sum\limits^{N-1}_{m=0} e^{j2\pi (k-k_0)m/N}
\end{align}$$

This orthogonal property assures that:
$$R(k) = \frac{1}{N} \sum\limits^{N-1}_{m=0} e^{j2\pi (k-k_0)m/N} = \delta(k_0)$$

Thus
$$z[k_0] = x[k_0] \times \sum\limits^{L-1}_{l=0} h_l e^{-j2\pi k_0 l/N} = x[k_0] \times H(k = k_0)\text{,}$$
where $H(k)$ is DFT of channel impulse reponse $(h_l, 0 \leq l \leq L-1)$ : is this the desired one-tap equalizer ?

Could someone tell me in which step I did use the condition flat fading $\Delta f < B_c \sim 1/\tau_m$ to come to single-tap equalizer model ?

Your mathematical derivation is correct, your $H[k]$ is the single-tap equalizer (i.e. one tap for each subcarrier, and the subcarriers do not mix with each other. That's the orthogonal in OFDM).

Let me try to explain this a bit more general, without going into coherence bandwidth and flat fading. To my understanding, explaining it with $B_c$ and flat fading is a bit superficial and undergraduate explanation, such that you understand the idea, but can't really prove it works. Even, the coherence bandwidth is not strictly defined, how can one design a system according to that?

Assuming that the channel is time-invariant over all times, the channel is an LTI system. The eigenfunctions of LTI systems are complex exponentials (of all frequencies). Note that these complex exponential range over all times (i.e. they are infinitely long). So, in principle, one could transmit data on all frequencies (infinitely small apart), but one would have to wait infinitely long until the signals arrived. The principle bases on the convolution theorem, i.e. your channel in time-domain does a convolution, so in frequency domain, it's doing elementwise multiplication. That's the continuous-time principle of OFDM. Clearly, not very attractive for practical implementations.

What can you do? We've know that two complex exponentials of different frequency are orthogonal to each other when considering the whole time axis, i.e. $\int_\mathbb{R}\exp(j2\pi f_1 t)\exp(j2\pi f_2 t)dt=\delta(f_1-f_2)$. We want to have orthogonality between the signals, (because we want to have single-tap equalization), but we dont want to wait infinitely long. Unfortunately, the orthogonality between to arbitrary frequency complex exponentials is not given, when the time interval is shorter: $\int_T\exp(j2\pi f_1 t)\exp(j2\pi f_2 t)dt\neq\delta(f_1-f_2)$. However, there are frequencies $f_1, f_2$ such that the orthogonality holds, namely $\int_T\exp(j2\pi f_1 t)\exp(j2\pi (f_1+\frac{n}{T})t)dt=\delta_n$ (here $\delta_n$ is the Kronecker symbol, i.e. 1 for $n=0$, else zero.). I.e. two exponentials that are $n/T$ apart in frequency are orthogonal over a time interval of $T$.

So, in principle, we could use time-intervals of $T$ and use only frequencies that are $1/T$ apart from each other to transmit orthogonal signals. But, there's more. If the orthogonality should hold at the receiver, the received signal needs to be a linear combination of complex exponentials. We know, that complex exponentials are eigenfunctions to the LTI system, so everything should be fine, right? No! Only infinitely long complex exponentials are eigenfunctions. If we transmit $\exp(j2\pi f_1t)\text{rect}(t/T)$, we will not get a complex exponential at the output, so it wont work. So, what to do? Essentially, we'd have to transmit infinitely long complex exponentials with frequency distance $1/T$ to get signals at the receiver, which are orthogonal in the interval $T$. Hm, still doesn't sound very appealing.

Here's the solution: The impulse response of the LTI system (i.e. the channel) is usually not infinitely long, but has a length $T_C$ (at least, we approximate/assume that it's zero for $t>T_C$). If the channel has length $T_C$, then for a given point in time $t$, only the signal times from $t-T_C$ to $t$ have an influence onto the output signal at the current point in time. So, we dont need to transmit infinitely long complex exponentials, but just of length $T+T_C$, and we get orthogonality in the time interval $[T_C,T+T_C]$.

You see where this is leading to, right? The extra time in the beginning is what we denote as the Cyclic prefix. Why do we call it cyclic? Because, the exponentials with frequency distance $1/T$ are periodic with period $T$. So, the signal in time $[0,...,T_C]$ is exactly equal to the signal at $[T,...,T+T_C]$.

To summarize: Transmitting complex exponentials of duration $T+T_C$ over an LTI system with impulse duration $T_C$ yields orthogonal signals over the interval $[T_C,...,T_C+T]$. Orthogonal means single-tap equalization is possible. What is the coefficient for the equalizer? Due to the convolution theorem, it is just the value of the frequency response of the channel at the carrier frequency.

This is the time-continuous explanation, which can become intuitive. However, I personally like more the discrete-time explanation, using linear algebra.

Let $F$ be the N-point Fourier transform matrix, i.e. $FF^H=F^HF=I$. So, the transmitted signal without CP is $y=F^Hx$ where $x$ is the data. Adding a CP is done by $y_{cp}=Cy$ with $C$ being a $N+N_{CP}\times N$ matrix. Then, we the received signal becomes

$$z=HCF^Hx$$ where $H$ is the Toeplitz matrix that does the convolution (I ignore noise here). $z$ is of length $N+N_{CP}+L-1$, where $L$ is the channel length. Then, at the receiver, the CP is removed, by a matrix $D$ of dimension $N\times N+N_{CP}+L-1$, i.e.

$$w=DHCF^Hx$$

Finally, we go to the frequency domain:

$$W=FDHCF^Hx$$

Let's have a look at the $N\times N$ matrix $DHC$: It's a matrix that performs circular convolution with the impulse response of your channel. In finite-discrete time the Convolution theorem states, that circular convolution equals elementwise multiplication of the DFTs. So, $FDHCF^H$ is a diagonal matrix, containing the frequency response of the channel on the diagonal. There you have your single-tap equalizer again (inversion by diagonal matrix is just division by the diagonal).


Can I simultaneously remove salt and pepper noise and improve the image resolution? What should my approach be towards designing a filter that can simultaneously remove salt & pepper noise and improve the resolution?

If you use two linear filters, they could be combined (because of linearity, the convolution with the convolution of the two filters equals to the consecutive convolution with the two filters).

$(I*a)*b=I*(a*b)$

Nevertheless the simplest way to deal with salt & pepper noise is median filtering, which is non-linear. If you use linear filters it is possible.


Check if the system is linear The system:
 $$ T(x[n]) = ax[n] + bx[n-3] $$

For me it seems that the system is linear:

$$
\begin{align}
T(\alpha_1x_1[n]+\alpha_2x_2[n]) & = a(\alpha_1x_1[n]+\alpha_2x_2[n]) + b(\alpha_1x_1[n-3]+\alpha_2x_2[n-3]) \\
& = \alpha_1(ax_1[n] + bx_1[n-3]) + \alpha_2(ax_2[n] + bx_2[n-3]) \\
& = \alpha_1T(x_1[n])+\alpha_1T(x_2[n])
\end{align}$$

Thus it's linear, however in the presentation I got it says it's not linear. (without reasoning) Where I'm making the mistake (or maybe there is a mistake in the presentation)

I believe there's either a mistake in the presentation or the presentation is using a different definition of linear.

For example, the system is linear in $x$ from a system perspective, but it's affine in $x[n]$ (and, therefore not linear) because of the $bx[n-3]$ offset.

On this site, we tend to go with the system definition rather than split hairs about linearity versus affine-ness.


Audio DSP HW Recommendation I'm very new to DSP hardware, so I am taking it upon myself to learn it with regard to Audio and Acoustic Signal Processing. I've learned quite a bit of theory in a graduate DSP course, but absolutely no practical implementation was introduced.

I've stumbled upon a site of miniDSP kits as well as a Texas Instruments TMS320C6748 Dev kit. What hardware should I purchase for learning audio DSP? And, if any, are there any good books that highlight the relative importance of different DSP hardware attributes? I would love to learn as much as I can about developing DSP hardware.

I think you might want to have a look at the stuff from XMOS. They have nice solutions (hardware, software, IDE) for developing audio applications on real(time) hardware. And besides that also nice development boards.


Do we need Channel estimation prior to (time domain DFE or linear) adaptative equalization in single tone transmission? If we are to transmit a single tone signal over a channel impaired by inter symbol interference and awgn, and if the equalization is adaptative (LMS/RLS track) and is done in time domain by LE (linear equalizer) or DFE, do we need to compute the channel impulse response first or the equalizer simply track channel changes when it computes it coefficients ?

Linear Equalizer: Depdends on the type, but no, generally you just estimate an optimal inverse ("optimal" according to the specific metric of linear equalizer type) of the channel impulse response and apply that.
Decision-Feedback Equalizer: No. The way the channel equalization is computed is literally in the name.


Time variance of $y(t)=x(2t)$ At first  I had trouble proving the time variance for this as I would do the usual mistake as in question Is the system represented by the equation $y(t) = x(2t)$ time invariant?. The solution has graphical representation of all the shifts in the process and they are the following.

I get the shift from the first picture to the second. It gets compressed, I know that property. I also understand what happens in the third picture. The $x_1(t)$signal is shifted to the right because of $t-2$ . The fifth picture is fine as well, it's just a shift to the right again from the second picture.

What I'm having trouble with is the 4th picture. I know that $y_2(t)=x_2(t)=x_1(2t-t_0)$ I forgot to mention, all the graphs are for $t_0=2$.

Shouldn't the 4th picture be a shift of 2 to the right?
And the last one maybe $(2t-2t_0)$ so $2t_0=4$ units to the right?

The signals plotted on the figure are correct.

As your provided link already shows: $$y_d(t) = \mathcal{T}\{x_d(t)\} = \mathcal{T}\{x(t-d)\} = \mathcal{T}\{x(t-2)\} = x(2t-2) \not= y(t-2) $$ where $$ ~ y(t)=\mathcal{T}\{x(t)\} = x(2t) \longrightarrow y(t-2) = x(2(t-2)) = x(2t-4) $$

Looking at the progression on the provided figure, $x_1(t)$ is a pulse of domain $~ [-2,2] ~$ then $y_1(t)=\mathcal{T}\{x_1(t)\} = x_1(2t)$ will be a pulse of domain $~ [-1,1] ~$ which is obtained by compressing $x(t)$ by $2$ about the origin $~t=0~$

$x_2(t)$ is obtained from $x_1(t)$, delaying it by $2$ units to the right, hence $x_2(t)$ has a domain of $~[0,4]$.
Therefore $y_2(t)=\mathcal{T}\{x_2(t)\} = x_2(2t)$ will be a pulse of domain $~ [0,2] ~$ which is obtained by compressing $x_2(t)$ by $2$ about the origin $~t=0~$

On the other hand shifting the pulse of $y_1(t)$ with the domain $[-1,1]$ by $2$ units to the right creates a new pulse with domain $[1,3]$ which is not equivalent to $[0,2]$; the domain of $y_2(t)$, from which we deduce, also graphically, that $y_1(t-2) \not = y_2(t)$


Causal system, order of numerator and denominator Is it true that if the order of numerator of a system function is larger than its denominator, then we can say that the system is causal? (Without finding R.O.C.)

Is there similar thing for Stability?

When $H(z)$ is rational, then the system is causal if and only if its ROC is the exterior of a circle outside the out-most pole, and the order of numerator is no greater than the order of the denominator.

A causal LTI system with a rational transfer function $H(z)$ is stable if and only if all poles of $H(z)$ are inside the unit circle of the z-plane, i.e., the magnitudes of all poles are smaller than 1.

so you also need to find the R.O.C.


Adaptive equalization vs inverse of transfer function I have the following equalization problem as shown in the figure below:

Now I can compute the coefficients for my adaptive FIR filter c (dim(c) = N) the following:
$\mathbf{c_{opt}} = (\mathbf{H}^T\mathbf{H})^{-1}\mathbf{H}~\mathbf{h_{ideal}}$
where $\mathbf{H}$ is a convolution matrix with shifted vectors of $\mathbf{h}$  and $\mathbf{h_{ideal}}$ is chosen such that $x[n]=d[n]$ (delay-free equalizer).

The channel impulse response is given as
$\mathbf{h} = [1, 0.5]^T$
$\Rightarrow H(z) = 1+0.5 z^{-1}$ so the inverse of the system would be IIR:
$1/H(z) = \frac{z}{z+0.5}$

Now the question is the following: What is the difference between the LS-solution with an adaptive filter and direct inversion of the system? Is it just that one filter is FIR and the other one IIR? Therefore with the FIR-filter we cannot reach full equalization and a residual error stays?

Inverting a channel can only be done when the channel is a minimum phase system (trailing echos only). A minimum phase system is characterized as having all zeros in the left half plane (for the s plane, or equivalently in a sampled system and the z plane all zeros inside the unit circle). Inverting such a channel results in poles where every zero exists, and a causal system that has any poles in the right half plane (outside the unit circle) is not stable.  So a minimum phase system has a stable causal inverse, while a mixed phase or maximum phase system does not.


trying to generate a sweep sine wave that increases its frequency till an upper limit and then starts decreasing its frequency to the lower limit Im trying to generate a sweep sine wave that increases its frequency till an upper limit and then starts decreasing its frequency to the lower limit at the same rate. Im trying to simulate in matlab.

The frequency of the wave will change based on the following formula: f = fmax-a*abs(T - t); a in the formula equals (fmax-fmin)/T
So the frequency should first increase and then decrease after it reaches time T at the same rate.

The matlab code is:

freq1 = 20;
freq2 = 200;
fs = 44100;
endTime = 0.1;
t = 1/fs:1/fs:(endTime*2);
a = (freq2-freq1)/endTime;
f = freq2 -   a*abs(endTime - t);
sinwave = sin(2*pi*f.*t);
plot(t,sinwave)

It seems to display a very different result when we plot the value obtained. After T the plot suddenly reverses, changes to original value of frequency and again increases in frequency I am not able to understand this behavior as the frequency seems to give values that increase and then decreases. Please help me understand why the plot for waves shows a different value and how to correct this. (Sorry not able to include plot due to web site restrictions)

I think that in your simulation, you only calculate the instantaneous frequency, but then fail to compute the phase for the sin() function correctly.

This is how you could correct this issue:

% parameters of linear up/down-sweep
freq1 = 20;
freq2 = 200;
fs = 44100;

% length of a sweep in one direction
endTime = 0.1;

% number of samples of a sweep in one sweep direction
N = endTime * fs;

% instantaneous frequency at each point in time ...
% ... first increasing for N points, then decreasing for N points
inst_f = [linspace(freq1, freq2, N) linspace(freq2, freq1, N)];

% since (in continuous time) instantaneous frequency is derivative of
% phase, we have to compute the "integral" to get the phase for sin().
phi = 2 * pi * cumsum(inst_f) / fs;
sweep = sin(phi);

% plotting
t = linspace(0, endTime * 2, length(s));
plot(t, sweep);

% wavwrite(sweep, fs, 32, '~/sweep_delete.wav');


Need help understanding 16QAM constellation diagram Is the identification of each point on the constellation diagram arbitrary as long as the sender and receiver know the mapping?

I have seen some diagrams which order in an outward spiral, while others go left to right. Does it matter?

Yes and no.

The mapping is arbitrary as long as the receiver correctly determines which constellation point a symbol is.  If the receiver makes a mistake, though, it is most likely going to pick a "neighbor" constellation point (i.e. a constellation point that is only one spot away).  It is highly unlikely that a correctly implemented receiver will pick a constellation point far away from the transmitted point.  Communication systems try to mitigate the impact of the "neighbor" mistakes by minimizing the Hamming distance between neighbors.  A map where all of the nearest neighbors have a Hamming distance of 1 is called a Gray map.  If you use a Gray map almost all of the symbol errors will only result in one bit error.


Improve contrast when image has a flood of white balancing I have an image taken a night time.

I apply _EqualizeHist to it. But the result is a high level of whiteness running through the image.

This is the image:

I had asked this question in a slightly different way and was given a link to improve the white balancing of the image.

This was the original question

This is that link:

enter link description here

I have got as far as this in my code but now I am stuck as to what to do next.  I would be grateful if anyone can suggest the next step?

This is my code:

            Image<Bgr, byte> frame = new Image<Bgr, byte>(file);
            frame._EqualizeHist();
            Image<Hsv, byte> hsv = new Image<Hsv, byte>(file);
            Image<Gray,byte> maskLight= new Image<Gray,byte>(hsv.Width,hsv.Height);
            Image<Gray, byte> maskDark = new Image<Gray, byte>(hsv.Width, hsv.Height);

            Emgu.CV.CvInvoke.cvThreshold(hsv[2], maskLight, 128, 255, THRESH.CV_THRESH_BINARY);
            Emgu.CV.CvInvoke.cvThreshold(hsv[2], maskDark, 0, 127, THRESH.CV_THRESH_BINARY);

            Bgr bgrLightest;
            MCvScalar scalarLightest;
            frame.AvgSdv(out bgrLightest, out scalarLightest, maskLight);

            Bgr bgrDarkest;
            MCvScalar scalarDarkest;
            frame.AvgSdv(out bgrDarkest, out scalarDarkest, maskDark);

UPDATE:
I have tried to implement the suggestion kindly goven by sansuiso.

This is the code:

 Image<Bgr, byte> rgb = new Image<Bgr, byte>(@"D:\789.jpg");
            pictureBox1.Image = rgb.ToBitmap();

            for (int x = 0; x < 360; x = x + 24)
            {
                for (int y = 0; y < 288; y = y + 24)
                {
                    rgb.ROI = new Rectangle(0, 0, 360, 288);
                    int xx = x;
                    int yy = y;
                    if (x > 0)
                    {
                        xx = xx - 12;
                    }
                    if (y > 0)
                    {
                        yy = yy - 12;
                    }
                    rgb.ROI = new Rectangle(xx, yy, 24, 24); //360*288
                    rgb._EqualizeHist();
                }
            }
            rgb.ROI = new Rectangle(0, 0, 360, 288);
            pictureBox2.Image = rgb.ToBitmap();

And this is the resulting image:

Obviously, I am not doing this right.  Can anyone see if my code is wrong please?

Histogram equalization is a global operation which can result in some areas being adversely effected at the cost of the rest of the image looking better. Consider using the Retinex algorithm. It works very well at improving both local and global contrast enhancement. There's a flurry of papers on this topic and the algorithm is somewhat straightforward to implement. The most challenging part is the selection of parameters.

Here's a list of papers and here's one that's very useful.


strange phenomenon in spectrogram of Windowed Fourier Transform [edited] I am trying to understand how the Windowed Fourier transform works, so I used Matlab to create an artificial signal, and plotted the Windowed Fourier transform using a gaussian window.

$x(t) = \left\{\begin{matrix}
 \sin(2\pi100t)& \text{for }0\leq{t}<\frac{1}{4}\\
 \sin(2\pi50t)& \text{for }\frac{1}{4}\leq{t}<\frac{1}{2}\\
 \sin(2\pi25t)& \text{for }\frac{1}{2}\leq{t}<\frac{3}{4}\\
 \sin(2\pi10t)& \text{for }\frac{3}{4}\leq{t}\leq1
\end{matrix}\right.$

and $f_s=500$, i.e. the data points are sampled at every 1/500 second.

This is what I have obtained for my plot.

However, for the last 2 intervals, $\frac{1}{2}\leq{t}<\frac{3}{4}$, and $\frac{3}{4}\leq{t}\leq1$, we can see some form of bifurcating phenomenon (vertically). Why is this so?

This happens because your window is too short. I don't have access to a plotting tool right now, but imagine for a second a slowly varying sinusoid that you chop up into pieces, and these pieces are shorter than the period of your sinusoid. If you take the Fourier transform of each of these pieces, some of your chunks will capture more energy of this sinusoid than others because they will contain crest or through of the function, while others will contain regions close to 0.

If you increase your window size, you'll be able to catch several cycles of your slow sinusoids inside each window, and the periodicity will go away. The tradeoff is that your windows will capture larger lengths of time, so boundaries between events will be more blurred.

If this still doesn't make sense, I'll try to include some plots later.


The blue curve does not seem accurate and is a little non-smooth. However, I will give a single-pole filter to approximate it.

It is not difficult to show that the step response of the filter $$H(z)=\frac{1-p}{1-pz^{-1}}$$
is in the form of

$$s[n]=(1-p^{n+1})u[n]$$

To find the step response, it suffices to shift the signal given in the question to the left and divide it by $2$. Hence, the output at zero is $0.25956786/2= 0.1298$. To find the parameter $p$, substitute it in the given equation for $s[n]$ at $n=0$:

$$0.1298=1-p\Rightarrow p=0.8702$$
So the filter we were looking for is $$H(z)=\frac{0.1298}{1-0.8702z^{-1}}$$
More clearly, $b=0.1298$ and  $a=[1,-0.8702]$.

Here is a comparison in Matlab:

It can be seen that the first $10$ outputs are very close to the actual output. However, after the output becomes non-smooth, there is a small error.
Alternatively, we could use the term at $n=22$ for our reference. In such case, $$1.95277071/2=1-p^{23}\Rightarrow p=0.8497.$$ This will lead to an accurate approximation at large $n$, but an error at small $n$.

Using the code

y = [0.25956786, 0.48537242, 0.68182373, 0.85275161, 1.00148237, 1.1309067, 1.24353671, 1.34155607, 1.42686391, 1.50111103, 1.62198067, 1.67093897, 1.71355379, 1.81104326, 1.83550966, 1.85680735, 1.87534702, 1.89148581, 1.90553486, 1.9177649, 1.92841125, 1.93767929, 1.95277071]
fun = @(P)2*(1-P.^(1:23))-y;
p = lsqnonlin(fun,.5);

$p=0.8644$ gives the least squares (LS) solution to this problem which is compared to the original output in the following plot ($\text{SSE}=0.02796$):

Note that, if the forward filter is $(b,a)$, then $(a,b)$ is the inverse filter which gives you the original input. If the error in the step response cannot be fixed, you can apply the LS approach on the overall output (taking into account the inverse filter) to get a better result.


problem finding 16 point DFT using two 8 point FFT (Divide and combine algorithm) MATLAB Write MATLAB code that determines and plot the N-point Discrete Fourier Transform of x[n] defined by the following equations:

x[n]=0.5*pi*n      n=0:16

Compute and plot 16-point DFT using two 8-point FFTs and combining techniques?

I will only give insight how to work on above problem but I won't give you the matlab code it can solved using two methods.
Method_1:- using DIT FFT algorithm hear inputs given in bit reversal order and output is multiplied with twiddle factor$w_n^k,where  k=0,1,2..,7$

So you use the same for your problem.
Method_2:- using DIF FFT algorithm hear the inputs are given in normal order at the output the butterfly's used are 2X2 & output is taken in bit reversal order like below


1D FFT algorithm elements explanation? Can someone explain to me what each element of the 1D FFT algorithm are (i.e. N,x,n,j,k) ? From here, http://zone.ni.com/reference/en-XX/help/371361H-01/lvanls/fft/#details

(Can't post pictures on stackoverflow just yet.... but its the 1D FFT description)

Consider a MATLAB/OCTAVE implementation of 1D-DFT/FFT sum:
$$X[k] = \sum_{n=0}^{N-1} {x[n]e^{-j \frac {2\pi}{N} k n} } $$
where $n=0,...,N-1$ and $k=0,...,N-1$

Those identifiers are used to denote the following:

$N$: is the FFT length as in $N$-point FFT.

$x[n]$: is the discrete-time signal whose DFT $X[k]$ is to be computed.

$X[k]$: is the complex-valued DFT of the signal $x[n]$

$n$: is the discrete-time index for the signal $x[n]$ in the range $n=0,...,N-1$.

$k$: is the frequency sample index of $X[k]$, in the range $k=0,...,N-1$.

$j$: is the imaginary unit $\sqrt {-1}$ in signal processing terminology.

Note if the signal $x[n]$'s length is shorter than the FFT length, then zero padding is applied on $x[n]$ to make its length $N$. Otherwise if it's longer than $N$, then the samples after $N-1$ are ignored according to MATLAB/OCTAVE implementation.


Is it correct to assume that $(-1)^n = \cos(\pi n)$ while computing the DTFT of $(-1)^n$? $\cos(\pi n)$ fluctuates between $-1$ and $1$ depending on the values of $n$, and it would be the same as computing it with an exponential but the problem is that I just get part of the right answer...

$$(-1)^n = \cos(n\pi)=(e^{j\pi})^n=(e^{-j\pi})^n=e^{j\pi n} = e^{-j\pi n}$$

All of these would work. However, applying a DFT is bit tricky since it's directly at the Nyquist frequency and it technically violates the sampling theorem. You can certainly do the math but interpreting the result is not straightforward. For example if you phase shift $\cos(n\pi)$ to $\sin(n\pi)$ the magnitude of the FFT drops to zero, although you only time-shifted the signal.


How to interpolate complex array (FT)? I do Fourier Transform for an image, then consider line(slice) in FFT image. When that slice isn't parallel neither x- not y-axis, coordinates are not integer.
E.g. slice can hold indices {(0,0), (0.25,1), (0.5,2), (0.75,3), (1,4)...}.
When using nearest-interpolation, quality of whatever I do is pure.

What techniques can you suggest to improve interpolation in Fourier grid?

The usual way to do this is to interpolate using surrounding points.  There is a simple 6-point mechanism that uses:

This question on the Mathematics SE site has some more details. Note the answer's correction to the equation in the question.